MATLAB Live Script to supplement "The effect of dopamine transporter blockade on optical self-stimulation: behavioral and computational evidence for parallel processing in brain reward circuitry"
This document was not intended to stand alone; It is assumed that anyone working through this script has already read the manuscript, with particular attention to the derivation of the reward-mountain model in the supporting-information file.
Preamble
In the current experiment, rats worked for optical stimulation of midbrain dopamine neurons. Time allocated to reward seeking was measured as a function of the strength and cost of the stimulation. We call the resulting three-dimensional data structure the "reward mountain." Dopamine neurotransmission was perturbed by administration of a dopamine-transporter blocker, and the consequent displacement of the reward mountain was determined.
The reward-mountain model was developed to account for data from experiments on rats working for rewarding electrical brain stimulation (electrical intracranial self-stimulation: eICSS). This model incorporates information obtained over many decades of research on the neural circuitry that intervenes between the tip of the stimulating electrode and the behavioral effects of the stimulation. Although this work has been dogged by uncertainty about the identity of the directly stimulated neurons, painstaking psychophysical experiments have provided extensive information about the directly activated neurons responsible for the rewarding effect, the operating principles of the neural circuitry in which they are embedded, and the reward-seeking behavior generated by this circuitry. Of particular importance to the current study are experiments characterizing the spatiotemporal integration of the electrically evoked reward signals and the growth of the rewarding effect as a function of stimulation strength.
In contrast to the case of eICSS, the identity of the directly stimulated substrate responsible for optical intracranial self-stimulation (oICSS) of midbrain dopamine neurons is known - the stimulation specifically and directly activates the dopaminergic neurons. However, only rudimentary information is available to date about the growth or spatiotemporal integration of the optically induced reward signal. One purpose of this document is to explore what the current experiment reveals about spatiotemporal integration and reward growth in the neural circuitry underlying oICSS. Another purpose is to explore the implications of the current work for understanding how eICSS and oICSS are related at the level of neural circuitry. We show by means of simulation that the results of the current experiment pose grave difficulties for an intuitively appealing account we and others proposed previously: that eICSS arises from the indirect activation of the directly stimulated dopamine neurons that give rise to oICSS. These difficulties challenge widely held notions about the organization of brain circuitry underlying reward.
Preliminaries
Setting the graph2files variable to true will cause graphs to be written to external files. These files are stored in a subdirectory of the current folder called 'Figures.' Setting this variable to false will speed execution.
The show_graphics variable determines whether graphics are stored within this Live Script. This variable must be set to false to enable efficient editing of this Live Script. When this Live Script has been executed with show_graphics set to true, subsequent editing and task switching will be slowed unacceptably. Set show_graphics to true only in preparation for saving updated html and/or pdf copies of this file. To prepare for editing after this script has been stored with show_graphics set to true, reset this variable to false, re-run the script, and save it.
Setting the tabs2files variable to true will store information about equations, figures, functions, and symbols in external Excel files. Setting this variable to false may speed execution slightly.
Setting the saveWS variable to true will cause the Matlab workspace to be stored in a .mat file in the folder from which this live script was run.
Be sure to set the default_dir to the directory you wish to use as the default on your system.
As its name implies, the 'Imported_Figures' folder contains images that this Live Script will import from external files. Those images should be stored in a sub-folder of the folder from which this Live Script will be run. The name of that folder is the argument of the set_impfigdir command below.
Blocks of redundant code are included in various sections. Although this lengthens the Live Script, it reduces the memory load due to the accumulation of workspace variables. Given that most variables are cleared at the boundaries between major sections, the blocks of redundant code enable debugging by means of "Run Section," "Run and Advance," etc. In a future revision, efficiency may be increased by aggregating the variables assigned in these blocks in structures that are saved and re-loaded.
t_start = tic; % start timer for entire script
tic; % start timer for the current section
global graphs2files show_graphics
graphs2files = true; % save graphics to external files?
tabs2files = true; % save tables of equation #s, figure #s, function #s, and symbol #s?
saveWS = true; % save final workspace?
show_graphics = true; % display graphics in this Live Script?
version = 11; % used in the names of stored files
% define and set the default directory
default_dir = '~/Work/Research/papers/Published/Opto_GBR/Simulations';
if ~exist(default_dir, 'dir')
disp(strcat({'Default directory '}, default_dir, {' does not exist.'}));
return
end
cd(default_dir);
% Define the directory that will contain figures generated here & create it if necessary.
global FigDir
FigDir = set_figdir('Figures', graphs2files); % Sub-directory of default directory
% Define the directory containing stored images to be imported.
ImpFigDir = set_impfigdir('Imported_figures'); % Sub-directory of default directory
% Define the directory that will house the tables & create if if necessary.
tabdir = set_tabdir('Tables', tabs2files); % Sub-directory of default directory
The following technical section implements preliminary steps required to set up the simulations.
In order to set up the simulations, some code must be executed to define basic functions and initialize variables.
This document was designed to run on any version of Matlab from R2018a onwards without requiring installation of external function files or scripts. All functions are either built-in (supplied by the Mathworks as part of the standard Matlab installation) or local (i.e., defined at the end of this Live-Script document). This should allow this script to execute on any standard installation of the supported Matlab versions.
The following code initializes the tables that store information about the formatted equations, figures and symbols used in this script.
%% Initialization
Initialize the tools that build the tables of equations, figures, functions, and symbols. Load the numbers and descriptions of the pre-defined functions (see below, "Tools for this live script" and "Building blocks for the functions included in the simulations."
global eqn_num eqn_tab fig_num fig_tab fun_num fun_tab sym_num sym_tab
init_all;
Load the numbers and descriptions of the pre-defined functions (see below, "Tools for this live script" and "Building blocks for the functions included in the simulations."
Enter symbols defined in the introductory paragraphs above.
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "eICSS", "electrical intracranial self-stimulation");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "oICSS", "optical intracranial self-stimulation");
keepVars = who; % Store variables to be retained
toc
The mountain model
tic;
The reward-mountain model was developed to account for operant performance as a function of the cost and strength of an experimenter-controlled reward, within the context of eICSS studies. The task performed to obtain the rewarding stimulation entails putting time on a clock by holding down a lever; a reward is delivered when the cumulative time the lever has been depressed reaches an experimenter-defined criterion (the "price" of the reward). The measure of operant performance employed is time allocation: the proportion of trial time that the rat devotes to procurement of the rewarding stimulation ("work"). The model is derived in the main body of the accompanying manuscript and then adapted for application to oICSS. Papers cited here are listed in the reference section of the main body of the manuscript.
In this section, we define the functional building blocks that will be used to simulate the output of the reward-mountain model.
We distinguish between the "shell" and "core" of the model. The shell consists of the variables that are observed (time allocation), manipulated (pulse frequency, price), and controlled (stimulation parameters held constant, physical work required to hold down the lever, affordances of the test environment for leisure activites, such as grooming, resting, and exploring). The core consists of the functions that compute the intensity of the reward produced by the stimulation train and combine this value with the opportunity and effort costs entailed in its procurement, thus generating what we call "payoffs." A single function based on the generalized matching law translates the payoffs generated in the core into the time-allocation values that are manifest in the shell.
Please see the main body of the text for definitions, explanations, and details.
Shell→core functions
The frequency-following function
Solomon et al. (2015) studied frequency-following fidelity in rats working for rewarding electrical stimulation of the medial forebrain bundle (MFB). They showed that the following function provides a good description of the relationship between the induced frequency of firing in the directly stimulated neurons and the pulse frequency:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "FreqFolElec", "Frequency of firing as a function of electrical pulse frequency");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "fF", "frequency-following function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Fpulse", "pulse frequency");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Fbend", "parameter determining sharpness of bend in the FreqFol function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Fro", "the pulse frequency in the center of the roll-off region");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Ffiring", "induced firing rate of the directly stimuled neurons");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Kf", "unit-translation constant for the frequency-following function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "FilterFun", "F, Fbend, Fro",...
"Frequency-following function");
The following graph relates the induced firing frequency to the electrical pulse frequency:
logFelecBend = 1.3222; % from Solomon et al., 2015
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587; % from Solomon et al., 2015
FelecRO = 10^logFelecRO;
Felec = logspace(0,3,121);
Fmfb = FilterFun(Felec, FelecBend, FelecRO); % Compute firing rate using the frequency roll-off function
FF_graph = plot_freqFoll(Felec,Fmfb,'mfb',2,1000,2,1000); % see "Functions that plot graphs and set attributes" below
if show_graphics
FF_graph.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "FreqFolMFB", ...
"Induced firing frequency in the directly stimulated neurons as a function of pulse frequency");
This graph shows that the firing of the directly stimulated neurons subserving eICSS of the medial forebrain bundle can follow the pulse frequency faithfully up to very high values (>= 350 pulses per second); the response of the neurons flattens abruptly as pulse frequency is increased further. This high-fidelity frequency following is consistent with the view that the directly stimulated cells subserving the rewarding effect are non-dopaminergic neurons with myelinated axons (Shizgal, 1997; Bielajew & Shizgal, 1986; Gallistel, Yeomans & Shizgal, 1981; Bielajew & Shizgal, 1982; Shizgal et al., 1980).
In addition to finding firing frequencies corresponding to a given pulse frequency, we will also need to do the reverse: finding the pulse frequency required to produce a given firing frequency. Thus, we define a back-solution of the frequency-following function:
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fFbacksolved", "Back-solution of the frequency-following function");
The neural-recruitment function
The 
function translates the pulse duration and current into the number of electrically excited neurons:
function translates the pulse duration and current into the number of electrically excited neurons:
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fN", "Neural-recruitment function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "c", "chronaxie of strength-duration function for pulses");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "d", "pulse duration");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "fN", "Neural-recruitment function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "I", "current");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Kns", "Neuron-recruitment function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Kis", "Current-distance constant");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "N", "Number of activated first-stage neurons");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "rhoI", "rheobase: threshold current to excite a first-stage neuron with a pulse of infinite duration");
Although we define this function here, we do not implement it or use it in the simulations. Instead, an output value is chosen that, together with a value for the rheobase of the strength-duration function for trains, generates location-parameter values within the range observed in past studies. This function is included in order to summarize the relationships on which existing tests of the counter model have been based (e.g., Simmons & Gallistel, 1994).
The burst-duration function
For completeness, we include a function that maps the duration of the pulse train into the duration of the evoked increase in the firing of the first-stage neurons. In practice, we assume that the two are equal, as will be case when frequency-following fidelity is high.
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fD", "Train-burst equation");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Dburst", "Duration of stimulation-induced increase in firing");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Dtrain", "train duration");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "fD", "Dtrain",...
"Train-burst function");
The subjective-probability function
This is another dummy function. As we review below, the subjective probability that a reward will be delivered upon satisfaction of the response requirement appears to equal the objective probability when the latter is 0.5 or higher. In the case of the present oICSS experiment, the objective probability of reward is always one.
n.b. A lower-case "p" is employed in the symbols.
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fp", "Subjective-probability equation");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "ProbObj", "Objective probability that reward will be delivered upon satisfaction of the response requirement");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "ProbSub", "Subjective probability that reward will be delivered upon satisfaction of the response requirement");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "ProbSubFun", "ProbObj",...
"Subjective-probability function");
The subjective-price (opportunity-cost) function
As described in Solomon et al., (2017), the subjective-price function is defined as
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "SubjPriceFun", "Subjective-price function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Pobj", "objective price (opportunity cost)");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Psub", "subjective price");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "PsubBend", "transition parameter of subjective-price function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "PsubMin", "minimum subjective price");
and plotted below.
PsubBend = 0.5;
PsubMin = 1.82;
Pobj = logspace(-1,2,120);
Psub = PsubFun(Pobj, PsubBend, PsubMin);
SP_graph = plot_PsubFun(Pobj,Psub, graphs2files, FigDir);
if show_graphics
SP_graph.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "PsubFun", ...
"The subjective-price function");
In addition to finding subjective prices corresponding to objective ones, we will also need to do the reverse. Thus, we define a back-solution of the subjective-price function:
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fPbacksolved", "Back-solution of the subjective-price function");
The subjective effort-cost function
The physical work entailed in holding down the lever is transformed into the subjective rate of exertion by the following function:
We do not attempt to model the form and parameters of this function. The dots over 
and 
signify that we define these quantities as rates.
and signify that we define these quantities as rates.[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fphiW", "subjective-effort-cost equation for work");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "fphi", "subjective effort-cost function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "dotPhiObj", "objective work rate required to hold down the lever");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "dotPhiSub", "subjective rate of exertion entailed in holding down the lever");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Kphi", "unit-conversion scalar for the effort cost of work");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "fphi", "dotPhiObj, Kec",...
"subjective-effort function");
The effect of the drug, if any, on the rate of subjective exertion is defined as:
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fphiWdrug", "drug-induced change in the subjective effort cost of work");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "KecW", "proportional drug-induced change in the subjective rate of exertion required to hold down the lever");
In the vehicle condition, 
assumes an implicit value of one.
assumes an implicit value of one.Activities such as grooming and exploring also entail performance of physical work. Thus, the subjective effort-cost function is also applied to these activities:
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fphiL", "subjective-effort-cost equation for leisure");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "dotPhiObjL", "average objective work rate required to perform leisure activities");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "dotPhiSubL", "average subjective rate of exertion entailed in performing leisure activities");
As in the case of the subjective rate of exertion entailed in work, we allow for drug-induced modulation of the subjective rate of exertion entailed in performance of leisure activities.
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fphiLdrug", "drug-induced change in the subjective effort cost of leisure activities");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "KecL", "proportional drug-induced change in the subjective rate of exertion required to perform leisure activities");
In the vehicle condition, 
assumes an implicit value of one.
assumes an implicit value of one.Core functions
The reward-growth function for brain stimulation reward (BSR)
The reward-growth function translates the aggregate rate of stimulation-induced firing into the intensity of the rewarding effect. Gallistel's team used operant matching to describe this function (Gallistel & Leon, 1991, Leon & Gallistel, 1992, Mark & Gallistel, 1993; Simmons & Gallistel, 1994). Shizgal (2003) proposed the following form for this function:
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fRbsr", "reward-growth equation for BSR");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "fRbsr", "reward-growth function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "FfiringHM", "firing frequency that produces a reward of half-maximal intensity");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "FpulseHM", "pulse frequency that produces a reward of half-maximal intensity");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "g", "reward-growth exponent");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Krg", "output scalar of the reward-growth function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Rbsr", "reward intensity produced by Ffiring");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "hatRbsr", "normalized reward intensity produced by a pulse frequency of Fpulse");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "fRbsr", "F, Fbend, Fhm, Fro, g, Krg",...
"Logistic reward-growth function");
It follows from the above that
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fRbsrNorm", "normalized reward-growth equation");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "fRbsrNorm", "F, Fbend, Fhm, Fro, g",...
"Normalized reward-growth function");
To accommodate the predicted rescaling of the input to the reward-growth function for oICSS by dopamine-transporter blockade, a scalar, 
, is added to the reward-growth equation in the drug condition of the experiment:
, is added to the reward-growth equation in the drug condition of the experiment:
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "RGfunDrug", "reward-growth equation for the drug condition");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Kda", "scalar representing the boost in dopamine release due to transporter blockade");
Thus,
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "FhmDrug", "FhmDrug as a function of FhmVehicle");
The following graphs plot the reward-growth function for three values of each parameter:
logFelecBend = 1.3222;
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587;
FelecRO = 10^logFelecRO;
numF = 121; % Number of pulse frequencies in Felec vector
Felec = logspace(0,3,numF);
FelecMat = repmat(logspace(0,3,numF),3,1); % add two rows in preparation for graphing
Fhm1 = 10^1.6; % The number of elements in the parameter vector must equal the # of rows in FelecMat
Fhm2 = 10^1.8;
Fhm3 = 10^2.0;
FhmVec = [Fhm1; Fhm2; Fhm3];
FhmMat = repmat(FhmVec,1,numF); % Store the Fhm values in a matrix of the same size as FelecMat
gElec1 = 5; % The number of elements in the parameter vector must equal the # of rows in FelecMat
gElec2 = 10;
gElec3 = 20;
gElecVec = [gElec1; gElec2; gElec3];
gElecMat = repmat(gElecVec,1,numF); % Store the gElec values in a matrix of the same size as FelecMat
Krg1 = 10^-0.2; % The number of elements in the paramer vector must equal the # of rows in FelecMat
Krg2 = 10^0;
Krg3 = 10^0.2;
KrgVec = [Krg1; Krg2; Krg3];
KrgMat = repmat(KrgVec,1,numF); % Store the Krg values in a matrix of the same size as FelecMat
RelecMat_Fhm = fRbsr(FelecMat, FelecBend, FhmMat, FelecRO, gElec1, Krg2);
RelecMat_gElect = fRbsr(FelecMat, FelecBend, Fhm2, FelecRO, gElecMat, Krg2);
RelecMat_Krg = fRbsr(FelecMat, FelecBend, Fhm2, FelecRO, gElec1, KrgMat);
TitleStrSemi = 'semi-log';
TitleStrLogLog = 'log-log';
pnam = "F_{hm}";
fnam = "Fhm";
% The data to be plotted must be in columns. Thus, FelecMat and RelecMat are transposed.
RG_Fhm_semilog = plot_RG(FelecMat',RelecMat_Fhm',pnam,FhmVec,fnam,TitleStrSemi,'lin');
RG_Fhm_loglog = plot_RG(FelecMat',RelecMat_Fhm',pnam,FhmVec,fnam,TitleStrLogLog,'log');
RG_dual_Fhm = dual_subplot(RG_Fhm_semilog, RG_Fhm_loglog, 'RG_Fhm_semilog_loglog',...
graphs2files,FigDir);
if show_graphics
RG_dual_Fhm.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "RGfunsElecFhm", ...
"Growth of reward intensity at three values of the position parameter");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "plot_RG", ...
"Fmat,RMat,pnam,pVec,fnam,TitleStr,linlog, varargins",...
"function to plot a single reward-growth function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "dual_subplot", ...
"g1, g2, dual_sub_out, graphs2files, figdir",...
"Function to plot two graphs side-by-side");
The semi-log plot is on the left, and the log-log plot is on the right.
These graphs shows how the position parameter rescales the input required to drive reward intensity to a particular level, thus sliding the reward-growth functions laterally.
pnam = "g";
fnam = "g";
RG_gElec_semilog = plot_RG(FelecMat',RelecMat_gElect',pnam,gElecVec,fnam,TitleStrSemi,'lin');
RG_gElec_loglog = plot_RG(FelecMat',RelecMat_gElect',pnam,gElecVec,fnam,TitleStrLogLog,'log');
RG_dual_g = dual_subplot(RG_gElec_semilog, RG_gElec_loglog, 'RG_gElec_semilog_loglog',...
graphs2files, FigDir);
if show_graphics
RG_dual_g.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "RGfunsElecG", ...
"Growth of reward intensity at three values of the steepness parameter");
The semi-log plot is on the left, and the log-log plot is on the right.
These graphs shows how the steepness parameter rotates the reward-growth functions around their midpoint.
pnam = "K_{rg}";
fnam = "Krg";
RG_Krg_semilog = plot_RG(FelecMat',RelecMat_Krg',pnam,KrgVec,fnam,TitleStrSemi,'lin');
RG_Krg_loglog = plot_RG(FelecMat',RelecMat_Krg',pnam,KrgVec,fnam,TitleStrLogLog,'log');
RG_dual_Krg = dual_subplot(RG_Krg_semilog, RG_Krg_loglog, 'RG_Krg_semilog_loglog',...
graphs2files, FigDir);
if show_graphics
RG_dual_Krg.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "RGfunsElecKrg", ...
"Growth of reward intensity at three values of the output-scaling parameter");
These graphs show how the parameter, 
, rescales the output of the reward-growth function vertically, thus determining its asymptotic level.
, rescales the output of the reward-growth function vertically, thus determining its asymptotic level. The preceding graphs show how rescaling of the input to and output from the reward-growth function produce orthogonal changes. This is true as long as as the firing of the directly stimulated neurons keeps pace with the pulse frequency. However, once the value of the position parameter nears the maximum firing frequency of the directly stimulated neurons (e.g., when the current is very low), rightward shifts of the reward-growth function, such as those produced by further decreases in current, are accompanied by decreases in its upper asymptote. The neurons can no longer fire fast enough to drive reward intensity to the same maximum as was achieved when the value of the position parameter was lower (e.g., because the current was higher). This has implications for the interpretation of changes in the location parameters of the reward mountain, as discussed below.
Fhm4 = 10^2.0;
Fhm5 = 10^2.3;
Fhm6 = 10^2.6
Fhm6 = 398.1072
FhmVec = [Fhm4; Fhm5; Fhm6];
FhmMat = repmat(FhmVec,1,numF); % Store the Fhm values in a matrix of the same size as FelecMat
RelecMat_FhmHi = fRbsr(FelecMat, FelecBend, FhmMat, FelecRO, gElec1, Krg2);
TitleStrSemi = 'HiF-semi-log';
TitleStrLogLog = 'HiF-log-log';
pnam = "F_{hm}";
fnam = "FhmHi";
% The data to be plotted must be in columns. Thus, FelecMat and RelecMat are transposed.
RG_FhmHi_semilog = plot_RG(FelecMat',RelecMat_FhmHi',pnam,FhmVec,fnam,TitleStrSemi,'lin');
RG_FhmHi_loglog = plot_RG(FelecMat',RelecMat_FhmHi',pnam,FhmVec,fnam,TitleStrLogLog,'log');
RG_dual_FhmHi = dual_subplot(RG_FhmHi_semilog, RG_FhmHi_loglog, 'RG_FhmHi_semilog_loglog',...
graphs2files,FigDir);
if show_graphics
RG_dual_FhmHi.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "RGfunsElecFhmHi", ...
"Growth of reward intensity at three values of the position parameter near FelecRO");
When the middle of the frequency roll-off zone (
) is positioned at 362 pulses 
(the median value reported by Solomon et al., 2015) and 
is 100 pulses , the normalized reward-growth function approaches an upper asymptote of one, as expected. However, doubling to ~200 pulses moderately decreases the upper asymptote. A further doubling to ~400 pulses pushes past , markedly truncating the growth of reward intensity.
) is positioned at 362 pulses (the median value reported by Solomon et al., 2015) and is 100 pulses , the normalized reward-growth function approaches an upper asymptote of one, as expected. However, doubling to ~200 pulses moderately decreases the upper asymptote. A further doubling to ~400 pulses pushes past , markedly truncating the growth of reward intensity.close all;
The location parameter of the reward-growth function for BSR
The location parameter of the reward-growth function is the firing rate that drives reward intensity to half its maximal value. The value of this parameter depends on the number of stimulated first-stage neurons, N, and the interval during which the stimulation train elevates their firing rate, 
.
. A prior study of temporal integration in the neural circuitry responsible for eICSS of the MFB (Sonnenschein et al., 2003) implies the following form for the function that determines 
:
:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "fHfiring", "strengh-duration function for trains expressed as a firing frequency");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "fHfiring", "C, D, N, RhoPi, varargin",...
"strengh-duration function for trains expressed as a firing frequency");
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "fH", "strengh-duration function for trains expressed as an aggregate firing frequency");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "fH", "C, D, RhoPi",...
"strengh-duration function for trains expressed as an aggreagate firing frequency");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Cmfb", "chronaxie of the strength-duration function for trains");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Dtrain", "duration of an electrical pulse train");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "RhoPi", "rheobase of the strength-duration function for trains");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "fH", "C, D, rhoPi",...
"Strength-duration function for trains; generates aggregate firing rate required to produce half-maximal reward intensity");
We assume that 
.
.The counter model is implicit in the equation for . According to this model, the aggregate rate of firing in the directly stimulated neurons determines the intensity of the rewarding effect.
. According to this model, the aggregate rate of firing in the directly stimulated neurons determines the intensity of the rewarding effect.The payoff from work
In keeping with the generalized matching law (Killeen, 1972), the benefit from work and its costs are combined in scalar fashion to yield a net payoff. We can format the expression for the payoff as a ratio of two rates (see: Gallistel & Gibbon, 2000):
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fUW", "payoff from work");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "fU", "utility function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "UW", "payoff from work");
or as a benefit/cost ratio:
The reward-growth function for leisure activities
We define a reward rate for the leisure activities that compete with pursuit of BSR as follows:
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fRaa", "Leisure reward rate");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Kaa", "subjective reward-rate scalar");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "DotRaa", "Average subjective reward rate from leisure activities");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "DotHatRaa", "normalized average subjective reward rate from leisure activities");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "fRaa", "Raa,Kaa",...
"reward-growth function for leisure activities");
The payoff from leisure activities
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, "fUL", "payoff from leisure");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "UL", "payoff from leisure activities");
The core → shell function
The behavioral-allocation function
The payoffs from pursuit of BSR and engagement in leisure activities are used by the sole core→shell function to compute the allocation of time to pursuit of BSR.
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "a", "payoff-sensitivity exponent");
The exponent 
determines how abruptly time allocation changes as a function of changes in payoff.
determines how abruptly time allocation changes as a function of changes in payoff.Restating this equation symbolically, we obtain
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "TU", "Time allocation defined in terms of payoffs");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "T", "time allocation");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Tmax", "maximum time allocation");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Tmin", "minimum time allocation");
We also define a normalized measure of time allocation
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "Tnorm", "Normalized time allocation");
as well as a value of time allocation, 
, at which the payoffs from work and leisure are equal, and thus, time allocation falls halfway between its minimal and maximal values:
, at which the payoffs from work and leisure are equal, and thus, time allocation falls halfway between its minimal and maximal values:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "Tmid", "Mid-range time allocation");
When the reward intensity produced by the stimulation approaches its upper asymptote, the price at which 
is:
is:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "Psub_e", "Definition of Psub_e");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "PsubEfun", ...
"dotPhiObj, Kaa, Kec, Krg, Raa, pObj, RnormMax, varargin",...
"function to compute PsubE");
To obtain an expression for 
, we back-solve the subjective-price equation, as shown above in
, we back-solve the subjective-price equation, as shown above indisp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='fPbacksolved')))}));
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "PobjEfun", ...
"dotPhiObj, Kaa, Kec, Krg, ObjAA, pObj, PsubBend, PsubMin, RnormMax, varargin",...
"Back-solution of the subjective-price function to return PobjE from PsubE");
We can now define the equation for normalized time allocation as:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "TnormExpand", "Normalized time allocation expanded");
and the full equation for time allocation as:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "Tfull", "Full 6-parameter time-allocation equation");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "TAfun", "a, F, Fbend, FpulseHM, Fro, g, Pobj, PobjE, PsubBend, PsubMin, RnormMax",...
"6-parameter TA function");
The reward mountain
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='Tfull')))}));
defines the surface of the reward-mountain model in a three-dimensional space defined by the two independent variables, the pulse frequency 
and the price 
, and a single dependent variable, time allocation
. Thus, this equation bridges the core of the model by linking three shell variables, two inputs and one output; the bridge is constructed from the shell → core, core, and core → shell functions. The slope of the mountain surface depends on the payoff-sensitivity parameter and the steepness parameter of the reward-growth equation, 
. The location of the mountain in the plane defined by the independent variables is determined by the two location parameters: whereas 
positions the mountain along the price axis, 
positions the mountain along the pulse-frequency axis. Finally, the altitude of the "valley floor" and "summit" are set by 
and 
, respectively.
and the price , and a single dependent variable, time allocation. Thus, this equation bridges the core of the model by linking three shell variables, two inputs and one output; the bridge is constructed from the shell → core, core, and core → shell functions. The slope of the mountain surface depends on the payoff-sensitivity parameter and the steepness parameter of the reward-growth equation, . The location of the mountain in the plane defined by the independent variables is determined by the two location parameters: whereas positions the mountain along the price axis, positions the mountain along the pulse-frequency axis. Finally, the altitude of the "valley floor" and "summit" are set by and , respectively.To simulate the mountain model, we must obtain values for the position parameters. . We proceed as follows:
. We proceed as follows:C = 0.473; % median from Sonnenschein et al., 2003
Dtrain = 0.5; % train duration in eICSS reward-mountain studies
logFelecBend = 1.3222; % from Solomon et al., 2015
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587; % from Solomon et al., 2015
FelecRO = 10^logFelecRO;
FPmax =1000; % for determining RbsrNormMax. This value is sufficiently above FelecRO to maximize the firing rate
N = 100; % arbitrary
RhoPi = 5000; % arbitrary
% RhoPi/N = 50, which is roughly consistent with Sonnenchein et al., 2003
FmfbHM = FpulseHMfun(C, Dtrain, FelecBend, FPmax, FelecRO, N, RhoPi)
FmfbHM = 97.3001
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "FpulseHMfun", ...
"C, D, Fbend, FPmax, Fro, N, RhoPi, varargin",...
"Function to compute the pulse frequency that produces half-maximal reward intensity");
The parameter that determines the location of the mountain along the price axis is .
. dotPhiObj = 1;
Kaa = 1;
Kec = 1;
Krg = 1;
dotRaa = 0.1;
pObj = 1;
PsubBend = 0.5;
PsubMin = 1.82;
Before computing , we must determine the maximum value of the normalized reward intensity (
). If the value of the position parameter of the reward-growth function is too high given the frequency roll-off parameters, 
will be less than one.
, we must determine the maximum value of the normalized reward intensity (). If the value of the position parameter of the reward-growth function is too high given the frequency roll-off parameters, will be less than one. [sym_num, sym_tab] = add_sym(sym_num, sym_tab, "RnormMax", "maximum normalized reward intensity");
g = 5;
%RnormMax = RnormMaxFun(Fmax, FelecBend, FmfbHM, FelecRO, g)
RnormMax = fRbsrNorm(FPmax, FelecBend, FmfbHM, FelecRO, g)
RnormMax = 0.9986
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RnormMax)
PobjE = 9.9860
We will now assign a value to the price-sensitivity parameter, a, and will simulate the mountain surface. (See Functions composing the reward-mountain model below.)
a = 3;
Felec = logspace(0,3,121)'; % column variable
Pobj = logspace(0,3,121); % row variable
Tmfb = TAfun(a, Felec, FelecBend, FmfbHM, FelecRO, g, Pobj, PobjE, PsubBend, PsubMin, RnormMax);
% n.b., numel(Felec) = numel(Pobj). Felec has been transposed. Thus, Tmfb is a square matrix.
MTN = plot_MTN(Felec, Pobj, Tmfb, 'off', 'MTN', 'reward mountain', ...
graphs2files, FigDir);
if show_graphics
MTN.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "Mountain", ...
"3D plot of the reward mountain");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "plot_MTN", ...
"F, Pobj, T, Visible, mtn_root, title_str, graphs2files, figdir, varargin",...
"Function to plot a single reward mountain");
clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
Moving the mountain: validation studies
tic;
close all;
We have carried out a series of experiments that test the mountain model. Given that there are two independent position parameters 
, it should be possible to move the mountain independently along either of the axes representing the objective price and the pulse frequency. This prediction has been tested in three ways.
, it should be possible to move the mountain independently along either of the axes representing the objective price and the pulse frequency. This prediction has been tested in three ways.Effect of varying the stimulation current
Changing the stimulation current alters the number of reward-related axons recruited. According to
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='FmfbHM')))}));
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='FelecHM')))}));
this will change the value of the parameter that positions the mountain along the pulse-frequency axis 
, but according to
, but according todisp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='Psub_e')))}));
and
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='fPbacksolved')))}));
changing the current will have no effect on the value of the parameter 
that positions the mountain along the price axis. We will now simulate the effect of changing the current by solving the time-allocation equation for two different values of N.
that positions the mountain along the price axis. We will now simulate the effect of changing the current by solving the time-allocation equation for two different values of N.C = 0.473; % median from Sonnenschein et al., 2003
D = 0.5; % typical train duration
logFelecBend = 1.3222; % from Solomon et al., 2015
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587; % from Solomon et al., 2015
FelecRO = 10^logFelecRO;
FPmax =1000; % for determining RbsrNormMax. This value is sufficiently above FelecRO to maximize the firing rate
N = [79,158]; % Number of neurons recruited by low & hi currents (2-element vector)
RhoPi = 5000; % arbitrary
FmfbHM = FpulseHMfun(C, D, FelecBend, FPmax, FelecRO, N, RhoPi) % FmfbHM is a 2-element vector
dotPhiObj = 1;
Kaa = 1;
Kec = 1;
Krg = 1;
dotRaa = 0.1;
pObj = 1;
PsubBend = 0.5;
PsubMin = 1.82;
FPmax = 1000; % This value is sufficiently above FelecRO to maximize the firing rate
g = 5;
RnormMax = fRbsrNorm(FPmax, FelecBend, FmfbHM, FelecRO, g) % RnormMax is a 2-element vector
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RnormMax)
% PobjE is a 2-element vector
a = 3;
Felec = logspace(0,3,121)'; % column variable
Pobj = logspace(0,3,121); % row variable
Tmfb1 = TAfun(a, Felec, FelecBend, FmfbHM(1), FelecRO, g, Pobj, PobjE(1), ...
PsubBend, PsubMin, RnormMax(1));
Tmfb2 = TAfun(a, Felec, FelecBend, FmfbHM(2), FelecRO, g, Pobj, PobjE(2), ...
PsubBend, PsubMin, RnormMax(2));
% n.b., numel(Felec) = numel(Pobj). Felec has been transposed. Thus, Tmfb is a square matrix.
MTNloI = plot_MTN(Felec, Pobj, Tmfb1, 'off', 'MTNloI', 'low current', ...
graphs2files, FigDir);
MTNhiI = plot_MTN(Felec, Pobj, Tmfb2, 'off', 'MTNhiI', 'high current', ...
graphs2files, FigDir);
dual_I_plot = dual_subplot(MTNloI, MTNhiI, 'MTNloI_hiI',...
graphs2files,FigDir);
if show_graphics
dual_I_plot.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "LoHiI_mtns", ...
"Effect of changing the current on the reward mountain");
(See Functions that plot graphs and set attributes below.)
Although the change in the position of the mountain following and increase in the number of reward-related neurons recruited can be discerned readily in the 3D plots, the effect is best depicted by means of contour graphs and a bar graph showing the changes in the position paramter. The contour graphs capture in two dimensions all of the spatial information in the 3D plots.
ContLoI = plot_contour(Felec, Pobj, Tmfb1, PobjE(1), FmfbHM(1), 'off', 'ContLoI', 'low current', ...
strcat({'N = '}, num2str(N(1))), graphs2files, FigDir);
ContHiI = plot_contour(Felec, Pobj, Tmfb2, PobjE(2), FmfbHM(2), 'off', 'ContHiI', 'high current', ...
strcat({'N = '}, num2str(N(2))), graphs2files, FigDir);
bg_LoHiI = plot_bg(FmfbHM(1), FmfbHM(2), PobjE(1), PobjE(2), 'off', 'bg_LoHiI',...
graphs2files, FigDir);
bg_root ='bg_LoHiI';
quad_I_plot = quad_subplot(ContLoI, ContHiI, bg_LoHiI, 'quad_LoHiI', bg_root, ...
graphs2files, FigDir);
if show_graphics
quad_I_plot.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "quad_LoHiI", ...
"Effect of changing the number of recruited neurons on the position of the reward mountain");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "plot_contour", ...
"F, Pobj, T, Pobj_e, Fhm, Visible, mtn_root, title_str, annot_str, graphs2files, figdir, varargin",...
"Function to plot the contour graph of a single mountain");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "quad_subplot", ...
"cont1, cont2, bg, quad_sub_out, bg_root, graphs2files, figdir",...
"Function to plot four graphs in a 2 x 2 mosaic");
The contour graph for the low-current condition is shown twice, once in the upper-left panel and once in the lower-right panel. By comparing the horizontal position of the mountain in the left column, the reader can quickly discern whether there has been a shift along the price axis, and by comparing the vertical position of the mountains in the bottom row, the reader can quickly discern whether there has been a shift along the pulse-frequency axis. As required by the mountain model, the latter shift is observed when the number of neurons recruited is increased due to boost in the stimulation current. The bar graph in the upper right provides a summary of the simulated shifts. The tiny rightward shift along the price axis is due to the fact that is a little closer to one at the higher current than at the lower current (due to the fact that is lower).
is a little closer to one at the higher current than at the lower current (due to the fact that is lower).Arvanitogiannis and Shizgal (2008) tested the effect of varying the stimulation current on the position of the reward mountain. In all four rats, the mountain shifted along the pulse-frequency axis as predicted. In two of these rats, there was no corresponding shift along the price axis, but rightward shifts were seen in the remaining two subjects. This deviation from the predictions was attributed to the heterogeneity of the stimulated neurons. Such a deviation would be expected if the stimulation activated two subpopulations of reward-related neurons that project to separate integrators with converging outputs (Arvanitogiannis, Waraczynski & Shizgal, 1996; Arvanitogiannis & Shizgal, 2008). The convergence model described below provides an example of such an arrangement.
clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
Effect of varying the train duration
tic;
close all;
In
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='quad_LoHiI')))}));
the mountain shifts along the pulse-frequency axis because changing the current alters the denominator of the ratio on the right-hand side of
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='fH')))}));
Changing the train duration should produce the same qualitative effect as changing the current, but by altering the numerator of the ratio on the right-hand side of
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='fH')))}));
instead. By increasing the duration of the train, there is more time for temporal summation and thus, the pulse frequency required to produce a reward of half-maximal intensity decreases.
C = 0.473; % median from Sonnenschein et al., 2003
D = 0.5; % typical train duration
logFelecBend = 1.3222; % from Solomon et al., 2015
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587; % from Solomon et al., 2015
FelecRO = 10^logFelecRO;
FPmax =1000; % for determining RbsrNormMax. This value is sufficiently above FelecRO to maximize the firing rate
N = 126;
D = [0.25,1.00]; % Short & long train durations (2-element vector)
RhoPi = 5000; % arbitrary
% RhoPi/N = 50, which is roughly consistent with Sonnenchein et al., 2003
FmfbHM = FpulseHMfun(C, D, FelecBend, FPmax, FelecRO, N, RhoPi) % FmfbHM is a 2-element vector
dotPhiObj = 1;
Kaa = 1;
Kec = 1;
Krg = 1;
dotRaa = 0.1;
pObj = 1;
PsubBend = 0.5;
PsubMin = 1.82;
g = 5;
RnormMax = fRbsrNorm(FPmax, FelecBend, FmfbHM, FelecRO, g) % RnormMax is a 2-element vector
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RnormMax)
% PobjE is a 2-element vector
a = 3;
Felec = logspace(0,3,121)'; % column variable
Pobj = logspace(0,3,121); % row variable
Tmfb1 = TAfun(a, Felec, FelecBend, FmfbHM(1), FelecRO, g, Pobj, ...
PobjE(1), PsubBend, PsubMin, RnormMax(1));
Tmfb2 = TAfun(a, Felec, FelecBend, FmfbHM(2), FelecRO, g, Pobj, ...
PobjE(2), PsubBend, PsubMin, RnormMax(2));
MTNloD = plot_MTN(Felec, Pobj, Tmfb1, 'off', 'MTNloD', 'short-duration train', ...
graphs2files, FigDir);
MTNhiD = plot_MTN(Felec, Pobj, Tmfb2, 'off', 'MTNhiD', 'long-duation train', ...
graphs2files, FigDir);
dual_D_plot = dual_subplot(MTNloD, MTNhiD, 'MTNLoHiD',...
graphs2files,FigDir);
if show_graphics
dual_D_plot.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "LoHiD_mtns", ...
"Effect of changing the train duration on the reward mountain");
ContLoD = plot_contour(Felec, Pobj, Tmfb1, PobjE(1), FmfbHM(1), 'off', 'ContLoD', 'short-duration train', ...
strcat({'D = '}, num2str(D(1))), graphs2files, FigDir);
ContHiD = plot_contour(Felec, Pobj, Tmfb2, PobjE(2), FmfbHM(2), 'off', 'ContHiD', 'long-duation train', ...
strcat({'D = '}, num2str(D(2))), graphs2files, FigDir);
bg_LoHiD = plot_bg(FmfbHM(1), FmfbHM(2), PobjE(1), PobjE(2), 'off', 'bg_LoHiD',...
graphs2files, FigDir);
bg_root ='bg_LoHiD';
quad_D_plot = quad_subplot(ContLoD, ContHiD, bg_LoHiD, 'quad_LoHiD', bg_root, ...
graphs2files, FigDir);
if show_graphics
quad_D_plot.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "LoHiD_quad", ...
"Effect of changing the train duration on the position of the reward mountain");
clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
tic;
close all;
Arvanitogiannis and Shizgal (2008) assessed the effect of increasing the train duration from 0.25 to 1.00 s in four rats. The data in all four cases correspond to the prediction: there were statistically reliable shifts of the reward mountain along the pulse-frequency axis but not along the price axis. The effect of increasing the train duration was tested in an additional six rats by Breton et al. (2014). In all six cases, the mountain shifted as predicted along the pulse-frequency axis, and in three of these cases, no reliable shifts along the price axis were observed. In the remaining three cases, increasing the train duration did produce rightward shifts along the price axis. These were again hypothesized to due to the recruitment of reward-related neurons projecting to separate integrators. The convergence model describe below offers a new interpretations of the rightward shifts.
if show_graphics
show_imported_graphic('TD_quad_Y14.png', 25, ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "LoHiD_quad_Y14", ...
"Rat Y14: Effect of changing the train duration on the position of the reward mountain");
A figure from Breton et al. (2014) showing single-subject data is reproducted above. Note the strong similarity between the simulated results in
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='LoHiD_quad')))}));
and the empirical results in
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='LoHiD_quad_Y14')))}));
toc
Effect of varying reward probability
tic;
close all;
Whereas changing the current or train duration is predicted to shift the reward mountain along the pulse-frequency axis but not along the price axis, changing the probability of delivering a reward upon satisfaction of the response requirement is predicted to produce an orthogonal shift: the mounain should move along the price axis but not along the pulse-frequency axis. This can be seen readily by inspection of
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='Pobj_e')))}));
,
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='StrDurFunTr')))}));
and
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='FmfbHM')))}));
. Intuitively, changing the reward probability should have no effect on the pulse-frequency required to produce a reward of half-maximal intensity, but is should rescale the payoff produced by this reward. Following the change in reward probability, the reward produced by a given pulse frequency is as intense as it was previously, but the ability of this pulse train to compete with alternate sources of reward will depend on the likelihood that the rat gets paid for the work it performs to obtain the electrical stimulation.
We will simulate the effect of changing the reward probability, first from 1.0 to 0.75 and then from 1.0 to 0.5:
C = 0.473; % median from Sonnenschein et al., 2003
D = 0.5; % typical train duration
logFelecBend = 1.3222; % from Solomon et al., 2015
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587; % from Solomon et al., 2015
FelecRO = 10^logFelecRO;
FPmax =1000; % for determining RbsrNormMax. This value is sufficiently above FelecRO to maximize the firing rate
N = 126;
D = 0.5; % train duration
RhoPi = 5000; % arbitrary
% RhoPi/N is roughly equal to 50, which is roughly consistent with Sonnenchein et al., 2003
FmfbHM = FpulseHMfun(C, D, FelecBend, FPmax, FelecRO, N, RhoPi)
FmfbHM = 77.2222
% n.b., FmfbHM1 = FmfbHM2
dotPhiObj = 1;
Kaa = 1;
Kec = 1;
Krg = 1;
dotRaa = 0.1;
pObj = [1, 0.75];
PsubBend = 0.5;
PsubMin = 1.82;
g = 5;
RnormMax = fRbsrNorm(FPmax, FelecBend, FmfbHM, FelecRO, g)
RnormMax = 0.9996
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RnormMax)
a = 3;
Felec = logspace(0,3,121)'; % column variable
Pobj = logspace(0,3,121); % row variable
Tmfb1 = TAfun(a, Felec, FelecBend, FmfbHM, FelecRO, g, Pobj, ...
PobjE(1), PsubBend, PsubMin, RnormMax);
Tmfb2 = TAfun(a, Felec, FelecBend, FmfbHM, FelecRO, g, Pobj, ...
PobjE(2), PsubBend, PsubMin, RnormMax);
MTNp1 = plot_MTN(Felec, Pobj, Tmfb1, 'off', 'MTNp1', 'p = 1.0', ...
graphs2files, FigDir);
MTNp0p75 = plot_MTN(Felec, Pobj, Tmfb2, 'off', 'MTNp0p75', 'p = 0.75', ...
graphs2files, FigDir);
dual_p0p75_plot = dual_subplot(MTNp1, MTNp0p75, 'MTNp1vsp0p75',...
graphs2files,FigDir);
if show_graphics
dual_p0p75_plot.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "p1vsp0p75_mtns", ...
"Effect of changing the reward probability on the reward mountain");
Contp1 = plot_contour(Felec, Pobj, Tmfb1, PobjE(1), FmfbHM, 'off', 'Contp1', 'p = 1.0', ...
strcat({'p = '}, num2str(pObj(1))), graphs2files, FigDir);
Contp0p75 = plot_contour(Felec, Pobj, Tmfb2, PobjE(2), FmfbHM, 'off', 'Contp0p75', 'p = 0.75', ...
strcat({'p = '}, num2str(pObj(2))), graphs2files, FigDir);
bg_p1vsp0p75 = plot_bg(FmfbHM, FmfbHM, PobjE(1), PobjE(2), 'off', 'p1vs0p75_bg',...
graphs2files, FigDir, -0.375, 0.075);
bg_root ='bg_p1vsp0p75';
quad_p0p75_plot = quad_subplot(Contp1, Contp0p75, bg_p1vsp0p75, 'quad_p1vsp0p75', bg_root, ...
graphs2files, FigDir);
if show_graphics
quad_p0p75_plot.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "quad_p1vsp0p75", ...
"Effect of changing the reward probability on the reward mountain");
clear -regexp ^bg ^Cont ^dual ^MTN ^quad; % Clear only restricted set of variables here
toc
tic;
close all;
pObj = [1, 0.5];
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RnormMax)
a = 3;
Felec = logspace(0,3,121)'; % column variable
Pobj = logspace(0,3,121); % row variable
Tmfb1 = TAfun(a, Felec, FelecBend, FmfbHM, FelecRO, g, Pobj, ...
PobjE(1), PsubBend, PsubMin, RnormMax);
Tmfb2 = TAfun(a, Felec, FelecBend, FmfbHM, FelecRO, g, Pobj, ...
PobjE(2), PsubBend, PsubMin, RnormMax);
MTNp1 = plot_MTN(Felec, Pobj, Tmfb1, 'off', 'MTNp1', 'p = 1.0', ...
graphs2files, FigDir);
MTNp0p5 = plot_MTN(Felec, Pobj, Tmfb2, 'off', 'MTNp0p5', 'p = 0.5', ...
graphs2files, FigDir);
dual_p0p5_plot = dual_subplot(MTNp1, MTNp0p5, 'MTNp1vsp0p5',...
graphs2files,FigDir);
if show_graphics
dual_p0p5_plot.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "p1vsp0p5_mtns", ...
"Effect of changing the reward probability on the reward mountain");
Contp1 = plot_contour(Felec, Pobj, Tmfb1, PobjE(1), FmfbHM, 'off', 'Contp1', 'p = 1.0', ...
strcat({'p = '}, num2str(pObj(1))), graphs2files, FigDir);
Contp0p5 = plot_contour(Felec, Pobj, Tmfb2, PobjE(2), FmfbHM, 'off', 'Contp0p5', 'p = 0.5', ...
strcat({'p = '}, num2str(pObj(2))), graphs2files, FigDir);
bg_p1vsp0p5 = plot_bg(FmfbHM, FmfbHM, PobjE(1), PobjE(2), 'off', 'p1vsp0p5_bg',...
graphs2files, FigDir, -0.375, 0.075);
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "p1vsp0p5_quad", ...
"Effect of changing the reward probability on the reward mountain");
bg_root ='bg_p1vsp0p5';
quad_p0p5_plot = quad_subplot(Contp1, Contp0p5, bg_p1vsp0p5, 'quad_p1vsp0p5', bg_root, ...
graphs2files, FigDir);
if show_graphics
quad_p0p5_plot.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "quad_p1vsp0p5", ...
"Effect of changing the reward probability on the reward mountain");
clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
tic;
close all;
Breton, Conover and Shizgal (2014) tested the effect of decreasing the reward probability from 1.00 to 0.75. The reward mountains obtained from all 10 rats shifted along the price axis as predicted. There was no statistically reliable shift along the pulse-frequency axis in the position of the mountains obtained from six of these rats. In the remaining four cases, the shifts along the pulse-frequency axis were small in comparison to the shifts along the price axis and were inconsistent in direction. Decreasing the probability of reward further to 0.5, increased the size of the shifts along the price axis in all seven rats tested. In five cases, there were no shifts observed along the pulse-frequency axis, and in the two cases in which such shifts were detected, they were inconsistent in direction. Data from one subject showing the effect of reducing the reward probability to 0.5 are shown below
if show_graphics
show_imported_graphic('prob_shifts_PD8.png',25,ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "p1vsp0p5_quad_PD8", ...
"Rat PD8: Effect of changing the reward probability on the reward mountain");
followed by a summary of the entire dataset.
if show_graphics
show_imported_graphic('prob_disc_all.png',30,ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "ProbDiscSummary", ...
"Summary of the effects of changing the reward probability on the reward mountain");
The upper panel of the summary shows that changing reward probability moves the mountain almost exclusively along the price axis and does so in a graded manner; the larger the change in reward probability, the larger the shift along the price axis. On the basis of
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='PayoffElec')))}));
we can estimate the subjective reward probability from the observed shift along the price axis. The lower panel of the summary shows that the subjective reward probabilities estimated in this manner are all but indistinguishable from the objective reward probabilities.
toc
The significance of orthogonal shifts
tic;
close all;
The validation experiments just described show conclusively that the reward mountain can be displaced along either the pulse-frequency or price axes, as predicted by the underlying model. In most cases, these displacements are orthogonal: the mountain shifts either along one axis or the other.
These observations have important implications for the form of the reward-growth function at the heart of the mountain model. In order for shifts along more than one axis to be possible, the form of the reward-growth function must distinguish changes that rescale its input and output. The logistic form implied by the matching data from the Gallistel lab instantiates such a distinction. Changes in the value of the position parameter, 
, rescale the input, thereby shifting the reward-growth function along the pulse-frequency axis, whereas changes in the value of the maximum attainable reward, , rescale the output, thereby shifting the maximum value of the reward-growth function vertically. In the figures below, we change the value of the input-rescaling position parameter by varying the current (thereby changing N, the number of reward-generating neurons recruited), and we rescale the output of the reward-growth function by changing the value of the maximum attainable reward, .
, rescale the input, thereby shifting the reward-growth function along the pulse-frequency axis, whereas changes in the value of the maximum attainable reward, , rescale the output, thereby shifting the maximum value of the reward-growth function vertically. In the figures below, we change the value of the input-rescaling position parameter by varying the current (thereby changing N, the number of reward-generating neurons recruited), and we rescale the output of the reward-growth function by changing the value of the maximum attainable reward, .logFelecBend = 1.3222; % from Solomon et al., 2015
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587; % from Solomon et al., 2015
FelecRO = 10^logFelecRO;
numF = 121; % Number of pulse frequencies in Felec vector
Felec = logspace(0,3,numF);
FelecMat = repmat(logspace(0,3,numF),3,1); % add two rows in preparation for graphing
C = median([0.453,0.503, 0.268, 0.224, 0.493, 0.64]); % from Sonnenschein et al., 2003
D = 0.5;
RhoPi = (5000 / fHfiring(C, D, 100, 5000)) * 100; % To position the logFhm values near 0.1 log10 units
gElec = 5;
NnLG1 = 10^1.8; % The number of elements in the parameter vector must equal the # of rows in FelecMat
NnLG2 = 10^2;
NnLG3 = 10^2.2;
NnLGvec = [NnLG1;NnLG2;NnLG3];
NnLGmat = repmat(NnLGvec,1,numF); % Store the N values in a matrix of the same size as FelecMat
KrgLG1 = 10^-0.2; % The number of elements in the parameter vector must equal the # of rows in FelecMat
KrgLG2 = 10^0;
KrgLG3 = 10^0.2;
KrgLGvec = [KrgLG1;KrgLG2;KrgLG3];
KrgLGmat = repmat(KrgLGvec,1,numF); % Store the Krg values in a matrix of the same size as FelecMat
FhmLG1 = fHfiring(C, D, NnLG1, RhoPi); % The number of elements in the parameter vector must equal the # of rows in FelecMat
FhmLG2 = fHfiring(C, D, NnLG2, RhoPi);
FhmLG3 = fHfiring(C, D, NnLG3, RhoPi);
FhmLGvec = [FhmLG1;FhmLG2;FhmLG3];
FhmLGmat = repmat(FhmLGvec,1,numF); % Store the Fhm values in a matrix of the same size as FelecMat
RelecLGmat_Nn = fRbsrFull(C,D,FelecMat,FelecBend,FelecRO,gElec,KrgLG2,NnLGmat,RhoPi);
RelecLGmat_Krg = fRbsrFull(C,D,FelecMat,FelecBend,FelecRO,gElec,KrgLGmat,NnLG2,RhoPi);
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "fRbsrFull", ...
"C, D, F, Fbend, Fro, g, Krg, N, RhoPi, varargin",...
"Full reward-growth function for BSR");
TitleStrSemiLog = 'logistic growth';
TitleStrLogLog = 'logistic growth';
pnam = "F_{hm}";
fnam = "Fhm";
% The data to be plotted must be in columns. Thus, FelecMat and RelecMat are transposed.
RG_NnLG_semilog = plot_RG(FelecMat',RelecLGmat_Nn',pnam,FhmLGvec,fnam,TitleStrSemiLog,'lin');
RG_NnLG_loglog = plot_RG(FelecMat',RelecLGmat_Nn',pnam,FhmLGvec,fnam,TitleStrSemiLog,'log');
pnam = "K_{rg}";
fnam = "Krg";
RG_KrgLG_semilog = plot_RG(FelecMat',RelecLGmat_Krg',pnam,KrgLGvec,fnam,TitleStrSemiLog,'lin');
RG_KrgLG_loglog = plot_RG(FelecMat',RelecLGmat_Krg',pnam,KrgLGvec,fnam,TitleStrSemiLog,'log');
dual_subplot(RG_NnLG_loglog, RG_KrgLG_loglog, 'RG_NnLG_KrgLG_loglog',...
graphs2files,FigDir);
if show_graphics
RG_NnLG_KrgLG_loglog.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "LogisticRGLL", ...
"When reward intensity grows as a logistic function of pulse frequency, shifts are orthogonal");
RG_NnLG_KrgLG_semilog = dual_subplot(RG_NnLG_semilog, RG_KrgLG_semilog, 'RG_NnLG_KrgLG_semilog',...
graphs2files,FigDir);
if show_graphics
RG_NnLG_KrgLG_semilog.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "LogisticRGSL", ...
"When reward intensity grows as a logistic function of pulse frequency, shifts are orthogonal");
The horizontal shifts of the reward-growth function are translated into shifts of the reward mountain along the pulse-frequency axis, whereas vertical shifts of the reward-growth function are translated in shifts of the reward mountain along the price axis. Above, we have illustrated such orthogonal shifts by simulating the effects of changing the current or train duration, as shown in
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='LoHiI_quad')))}));
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='LoHiD_quad')))}));
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='LoHiD_quad_Y14')))}));
and the effects of changing the reward probability, as shown in
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='p1vsp0p75_quad')))}));
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='p1vsp0p5_quad')))}));
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='p1vsp0p5_quad_PD8')))}));
What would happen if frequency-following fidelity remained the same, but the input-scaling and output-scaling parameters of the reward-growth function were no longer independent? We can simulate such a case by replacing the logistic reward-growth function with a power function:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "RGpg", "Power growth of reward intensity");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "fRpg", "F, Fbend, Fro, g, Kin, Krg",...
"Function to compute power growth of reward");
This power-growth function can be rewritten as
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "RGpgReformat", "Power growth of reward intensity");
In contrast to the case of logistic growth, the two scaling constants act jointly, in an inseparable manner, and exclusively on the output of the function. Changing either scaling constant shifts the power-growth function along the y axis but does not change its position along the x axis. As in the case of the logistic reward-growth function, g determines the steepness of reward growth (the slope on double-logarithmic coordinates). Although 
appears to scale the output, whereas
appears to scale the induced firing frequency,
appears to scale the output, whereas appears to scale the induced firing frequency, disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='RGpgReformat')))}));
shows that the two constants function as one. Thus, reward intensity will grow as a power function of the electrical pulse frequency until the induced firing frequency in the stimulation reward-related neurons asymptotes. Changing the current displaces the logistic reward-growth function laterally, a shift that is orthogonal to the one produced by changing the scale parameter, 
. In contrast, changing the current displaces the power reward-growth function vertically, the same direction as the shift produced by changing the scale parameter. This is shown graphically below.
. In contrast, changing the current displaces the power reward-growth function vertically, the same direction as the shift produced by changing the scale parameter. This is shown graphically below. In the following simulations, changing the value of Nserves as a proxy for changes in current. The parameter is inversely proportional to the number of stimulated neurons (N). In the figure legends, the value of the parameter is normalized to the value corresponding to the middle value of N.
parameter is inversely proportional to the number of stimulated neurons (N). In the figure legends, the value of the parameter is normalized to the value corresponding to the middle value of N.logFelecBend = 1.3222;
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587;
FelecRO = 10^logFelecRO;
FPmax = 1000;
numF = 121; % Number of pulse frequencies in Felec vector
Felec = logspace(0,3,numF);
FelecMat = repmat(logspace(0,3,numF),3,1); % add two rows in preparation for graphing
gElec = 2.5;
NnPG1 = 10^(2-(0.3/gElec)); % The number of elements in the parameter vector must equal the # of rows in FelecMat
NnPG2 = 10^2;
NnPG3 = 10^(2+(0.3/gElec));
% NnPGvec = [NnPG1;NnPG2;NnPG3];
% NnPGmat = repmat(NnPGvec,1,numF); % Store the N values in a matrix of the same size as FelecMat
FFaggNorm = FilterFun(FPmax,FelecBend, FelecRO) * NnPG2; % Normalization factor
% When Kout = 1 and NnPG = NnPG2, FFaggNorm ensures that R = 1 when FF = FFmax
KinPGnorm = FFaggNorm/NnPG2;
KinPG1 = FFaggNorm/NnPG1;
KinPG2 = FFaggNorm/NnPG2;
KinPG3 = FFaggNorm/NnPG3;
KinPGvec = [KinPG1;KinPG2;KinPG3];
KinPGmat = repmat(KinPGvec,1,numF); % Store the Kin values in a matrix of the same size as FelecMat
% The number of elements in the parameter vector must equal the # of rows in FelecMat
KoutPG1 = 10^-0.3;
KoutPG2 = 10^0;
KoutPG3 = 10^0.3; % RpgMax = 10^0.15 = 1.4125
KoutPGvec = [KoutPG1;KoutPG2;KoutPG3];
KoutPGmat = repmat(KoutPGvec,1,numF); % Store the Kout values in a matrix of the same size as FelecMat
RelecPGmat_Kin = fRpg(FelecMat,FelecBend,FelecRO,gElec,KinPGmat,KoutPG2);
RelecPGmat_Kout = fRpg(FelecMat,FelecBend,FelecRO,gElec,KinPG2,KoutPGmat);
TitleStrSemiLog = 'power growth';
TitleStrLogLog = 'power growth';
pnam = "K_{in} / K_{in_{norm}}";
fnam = "Kin";
KinLgndVec = KinPGvec ./ KinPGnorm;
RG_KinPG_semilog = plot_RG(FelecMat',RelecPGmat_Kin',pnam,KinLgndVec,fnam,TitleStrSemiLog,'lin');
RG_KinPG_loglog = plot_RG(FelecMat',RelecPGmat_Kin',pnam,KinLgndVec,fnam,TitleStrLogLog,'log');
pnam = "K_{out}";
fnam = "Kout";
RG_KoutPG_semilog = plot_RG(FelecMat',RelecPGmat_Kout',pnam,KoutPGvec,fnam,TitleStrSemiLog,'lin');
RG_KoutPG_loglog = plot_RG(FelecMat',RelecPGmat_Kout',pnam,KoutPGvec,fnam,TitleStrLogLog,'log');
RG_KinPG_KoutPG_semilog = dual_subplot(RG_KinPG_semilog, RG_KoutPG_semilog, 'RG_KinPG_KoutPG_semilog',...
graphs2files,FigDir);
if show_graphics
RG_KinPG_KoutPG_semilog.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "PowerRGSL", ...
"When reward intensity grows as a power function, all shifts are vertical");
dual_subplot(RG_KinPG_loglog, RG_KoutPG_loglog, 'RG_KinPG_KoutPG_loglog',...
graphs2files,FigDir);
if show_graphics
RG_KinPG_KoutPG_loglog.Visible = 'on';
end
%shg
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "PowerRGLL", ...
"When reward intensity grows as a power function, all shifts are vertical");
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='RGpg')))}));
Unlike the case with logistic growth of reward intensity, the variable that scales the input to the power reward-growth function 
and the variable that scales its output 
are no longer independent and cannot move the reward-growth curve in orthogonal directions. As a result, the reward mountain can shift only along a single axis when the underlying reward-growth function lacks independent input- and output-scaling parameters. This is demonstated in the following section.
and the variable that scales its output are no longer independent and cannot move the reward-growth curve in orthogonal directions. As a result, the reward mountain can shift only along a single axis when the underlying reward-growth function lacks independent input- and output-scaling parameters. This is demonstated in the following section.clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
tic;
close all;
The reward mountain can shift only along the price axis when the underlying reward-growth function lacks indepedent input- and output-scaling parameters
dotPhiObj = 1;
Kaa = 1;
Kec = 1;
dotRaa = 0.1;
pObj = 1;
a = 3;
Felec = logspace(0,3,121)'; % column variable
Pobj = logspace(0,3,121); % row variable
logFelecBend = 1.3222;
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587;
FelecRO = 10^logFelecRO;
gElec = 2.5;
N = 100;
PsubBend = 0.5;
PsubMin = 1.82;
FPmax = 1000;
FFmax = FilterFun(FPmax,FelecBend, FelecRO);
Kin = [FFmax, FFmax / 2^(1/gElec)];
Kout = [1,2];
RmaxBase = fRpg(1000,FelecBend,FelecRO,gElec,Kin(1),Kout(1));
RmaxKin = fRpg(1000,FelecBend,FelecRO,gElec,Kin,Kout(1));
RmaxKout = fRpg(1000,FelecBend,FelecRO,gElec,Kin(1),Kout);
PsubEpgBase = PsubEpgFun(dotPhiObj, Kaa, Kec, Kout, dotRaa, pObj, RmaxBase);
PsubEpgKin = PsubEpgFun(dotPhiObj, Kaa, Kec, Kout, dotRaa, pObj, RmaxKin);
PsubEpgKout = PsubEpgFun(dotPhiObj, Kaa, Kec, Kout, dotRaa, pObj, RmaxKout);
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "PsubEpgFun", ...
"dotPhiObj, Kaa, Kec, Krg, ObjAA, pObj, Rmax",...
"Function to compute PsubE for power reward growth");
TpgBase = TApgFun(a, Felec, FelecBend, FelecRO, gElec, Kin(1), Kout(1), Pobj, PsubBend, PsubEpgBase, PsubMin);
TpgKin2 = TApgFun(a, Felec, FelecBend, FelecRO, gElec, Kin(2), Kout(1), Pobj, PsubBend, PsubEpgKin(2), PsubMin);
TpgKout2 = TApgFun(a, Felec, FelecBend, FelecRO, gElec, Kin(1), Kout(2), Pobj, PsubBend, PsubEpgKout(2), PsubMin);
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "TApgFun", ...
"a, Felec, FelecBend, FelecRO, gElec, Kin, Kout, N, Pobj, PsubBend, PsubEpgBase, PsubMin",...
"Time allocation in response to power growth of reward intensity");
xmin = 0;
xmax = 2.3;
ymin = 1.5;
ymax = 3;
MTNpgKin1 = plot_MTN(Felec, Pobj, TpgBase, 'off', 'MTNpgKin1', 'Kin = Kin1', ...
graphs2files, FigDir, xmin, xmax, ymin, ymax);
MTNpgKin2 = plot_MTN(Felec, Pobj, TpgKin2, 'off', 'MTNpgKin2', strcat('Kin1','\div','2.0^{1/g}'), ...
graphs2files, FigDir, xmin, xmax, ymin, ymax);
dual_pgKin_plot = dual_subplot(MTNpgKin1, MTNpgKin2, 'MTNpgKin1Kin2',...
graphs2files,FigDir);
if show_graphics
dual_pgKin_plot.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "PG_K1_1vsK1_2_mtns", ...
"Effect of changing the power-growth input-scaling parameter");
MTNpgKout1 = plot_MTN(Felec, Pobj, TpgBase, 'off', 'MTNpgKout1', 'Kout = Kout1', ...
graphs2files, FigDir, xmin, xmax, ymin, ymax);
MTNpgKout2 = plot_MTN(Felec, Pobj, TpgKout2, 'off', 'MTNpgKout2', 'Kout = 2*Kout1', ...
graphs2files, FigDir, xmin, xmax, ymin, ymax);
dual_pgKout_plot = dual_subplot(MTNpgKout1, MTNpgKout2, 'MTNpgKout1Kout2',...
graphs2files,FigDir);
if show_graphics
dual_pgKout_plot.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "PG_Kout1vsKout2_mtns", ...
"Effect of changing the power-growth output-scaling parameter");
FhmPGKin = find_FhmPG(FPmax, FelecBend, FelecRO, gElec, Kin, Kout(1), RmaxKin);
FhmPGKout = find_FhmPG(FPmax, FelecBend, FelecRO, gElec, Kin(1), Kout, RmaxKout);
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "find_FhmPG", ...
"FPmax, Fbend, Fro, g, Kin, Kout, RmaxPG",...
"Function to compute the Fhm value for power reward growth");
PobjEpgBase = PsubBsFun(PsubEpgBase, PsubBend, PsubMin);
PobjEpgKin = PsubBsFun(PsubEpgKin, PsubBend, PsubMin);
PobjEpgKout = PsubBsFun(PsubEpgKout, PsubBend, PsubMin);
ContKin1 = plot_contour(Felec, Pobj, TpgBase, PobjEpgKin(1), FhmPGKin(1), 'off', 'ContKin1', 'Kin = 1.0', ...
strcat({'Kin = '}, num2str(round(Kin(1),2))), graphs2files, FigDir,...
xmin, xmax, ymin, ymax);
ContKin2 = plot_contour(Felec, Pobj, TpgKin2, PobjEpgKin(2), FhmPGKin(2), 'off', 'ContKin2', 'Kin = 2.0^{1/g}', ...
strcat({'Kin = '}, num2str(round(Kin(2),2))), graphs2files, FigDir,...
xmin, xmax, ymin, ymax);
bg_Kin1vskKin2 = plot_bg(FhmPGKin(1), FhmPGKin(2), PobjEpgKin(1), PobjEpgKin(2), 'off', 'Kin1vsKin2_bg',...
graphs2files, FigDir);
bg_root = 'bg_Kin1vsKin2';
quad_Kin1vsKin2 = quad_subplot(ContKin1, ContKin2, bg_Kin1vskKin2, 'quad_Kin1vsKin2', bg_root, ...
graphs2files, FigDir);
if show_graphics
quad_Kin1vsKin2.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "quad_Kin1vsKin2", "Effect of increasing the value of the input-scaling paramater");
ContKout1 = plot_contour(Felec, Pobj, TpgBase, PobjEpgKout(1), FhmPGKout(1), 'off', 'ContKout1', 'Kout = Kout1', ...
strcat({'Kout = '}, {'Kout1'}), graphs2files, FigDir,...
xmin, xmax, ymin, ymax);
ContKout2 = plot_contour(Felec, Pobj, TpgKout2, PobjEpgKout(2), FhmPGKout(2), 'off', 'ContKout2', 'Kout = 2 * Kout1', ...
strcat({'Kout = '}, {'2 * Kout1'}), graphs2files, FigDir,...
xmin, xmax, ymin, ymax);
bg_Kout1vskKout2 = plot_bg(FhmPGKout(1), FhmPGKout(2), PobjEpgKout(1), PobjEpgKout(2), 'off', 'Kout1vsKout2_bg',...
graphs2files, FigDir);
bg_root = 'bg_Kout1vsKout2';
quad_Kout1vsKout2 = quad_subplot(ContKout1, ContKout2, bg_Kout1vskKout2, 'quad_Kout1vsKout2', bg_root, ...
graphs2files, FigDir);
if show_graphics
quad_Kout1vsKout2.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "quad_Kout1vsKout2", "Effect of increasing the value of the output-scaling paramater");
The graphs shown above demonstrate that changes to the input- and output-scaling parameters of the power-growth function produce identical shifts in the position of the reward mountain. Changing the value of either parameter shifts the mountain rightward along the price axis, in sharp contrast to the orthogonal shifts observed when reward intensity grows as a logistic function of pulse frequency.
In the case of eICSS, we know from the experiments of Gallistel's group that reward intensity grows in a manner similar to a logistic function of pulse frequency. Unlike the case of power growth, logistic growth is characterized by independent parameters that scale the input and output. This characteristic is what causes the effect of varying reward probability to shift the reward mountain in a direction orthogonal to the effects of varying current or train duration. It follows that when such orthogonal shifts are observed, this requires that the underlying reward-growth function incorporates independent parameters that scale the input (the position parameter) and the output (the parameter representing the maximum reward intensity). Thus even in the absence of data from matching experiments or other methods for measuring the growth of reward intensity, we can make inferences about the reward-growth function by observing the way in which various manipulations shift the reward mountain. This will be important for the interpretation of the oICSS data.
clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
At what stage of processing does perturbation of dopaminergic neurotransmission alter reward seeking in eICSS?
tic;
close all;
We have described changes in the position parameter as alterations in the sensitivity of the reward-growth function. One might best label as the "inverse-sensitivity" parameter because the lower its value, the weaker the input required to drive reward intensity to a given level. In contrast, we have described changes in the output-scaling parameter as alterations in the gain of the reward-growth function. Just as turning the volume knob on an audio amplifier changes the loudness produced by weak and strong inputs by the same percentage, changing alters the reward intensity produced by both low and high pulse frequencies (within the frequency-following range) by the same percentage. Changes in sensitivity shift the reward mountain along the pulse-frequency axis, whereas changes in gain shift the mountain along the price axis.
as alterations in the sensitivity of the reward-growth function. One might best label as the "inverse-sensitivity" parameter because the lower its value, the weaker the input required to drive reward intensity to a given level. In contrast, we have described changes in the output-scaling parameter as alterations in the gain of the reward-growth function. Just as turning the volume knob on an audio amplifier changes the loudness produced by weak and strong inputs by the same percentage, changing alters the reward intensity produced by both low and high pulse frequencies (within the frequency-following range) by the same percentage. Changes in sensitivity shift the reward mountain along the pulse-frequency axis, whereas changes in gain shift the mountain along the price axis.Operant peformance in the eICSS paradigm is altered by perturbation of dopaminergic neurotransmission. Boosting dopaminergic signalling typically enhances performance, whereas attenuating dopaminergic signalling typically attenuates performance. These changes have long been attributed to the modulation of reward sensitivity by dopaminergic agents (See: Hernandez et al., 2010). If so, drugs that alter dopaminergic neurotransmission should shift the reward mountain along the pulse-frequency axis and not along the price axis. We have shown that this is not so. In rats working for electrical stimulation of the medial forebrain bundle, we demonstrated that enhancement of dopaminergic signalling by the highly specific reuptake blocker, GBR-12909, shifted the reward mountain rightward along the price axis in 8/10 rats without producing any statistically reliable shifts along the pulse-frequency axis (Hernandez et al., 2012). The D2/D4/5HT7 receptor blocker, pimozide, shifted the reward mountain leftward along the price axis in 5/6 rats without producing any statistically reliable shifts along the pulse-frequency axis (Trujillo-Pisanty, Conover & Shizgal, 2014). To account for these complementary effects, the changes in dopaminergic neurotransmission had to have altered one or more of the variables that determine the position of the reward mountain along the price axis, such as the output-scaling parameter of the reward-growth function , the subjective effort cost (), or the payoff from alternate activities (
).
, the subjective effort cost (), or the payoff from alternate activities ().if show_graphics
show_imported_graphic('GBR_pimozide_summary.png', 100, ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "PimGBRshifts", ...
"Changes in dopaminergic neurotransmission shift the mountain along the price axis");
An intuitively appealing interpretation of these data is based on the notion that the neurons subserving eICSS of the medial forebrain bundle project to midbrain dopamine neurons and that it is the transynaptic activation of the dopamine neurons that gives rise to the reward effect of the medial forebrain bundle stimulation (the "series" model of brain-reward circuitry). In the following section, we apply the reward mountain model to the oICSS data from the present paper. We show that the series model cannot account for the data from both the eICSS and oICSS studies. After demonstrating this, we discuss "convergence" models that show promise for an intergrated account of both eICSS and oICSS, and we propose experiments that put such accounts to empirical test.
toc
Optical intracranial self-stimulation of midbrain dopamine neurons
tic;
close all;
In the experiments reviewed above, rats worked to trigger trains of electrical current pulses applied to the medial forebrain bundle. In contrast, in the experiment we now discuss, rats worked to trigger trains of optical pulses delivered through an optical fiber positioned over the ventral tegmental area of the midbrain. As a result of viral transcription and Cre-Lox recombination, dopamine neurons with somata in, or adjacent to, the VTA expressed channelrhodopsin-2 (ChR2) and thus could be excited by light delivered at an appropriate wavelength (473 nm). As others have reported previously (Witten et al., 2012), rats expressing ChR2 in midbrain dopamine neurons work vigorously to obtain such optical stimulation.
We show that TH-Cre rats expressing ChR2 in midbrain dopamine neurons learn to perform the cumulative hold-down task to receive optical stimulation of the ventral tegmental area (VTA), where the somata of dopamine neurons projecting to forebrain targets reside. As in the case of eICSS, time allocation varies smoothly as a function of the strength and price of the optical stimulation. The surface defined by the mountain model fits the data well.
if show_graphics
show_imported_graphic('Bechr29_veh_mtn.png', 30, ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "BeChR29vehMtn", ...
"The mountain model fitted to vehicle-condition data from rat BeChR29");
The effect of the dopamine-transporter blocker, GBR-12909.
To boost dopaminergic neurotransmission, we administered 20 mg/kg of the dopamine-transporter blocker, GBR-12909. Under the influence of the drug, the mountain fitted to the data from rat BeChR29 shifted leftwards along the pulse-frequency axis and rightwards along the price axis.
if show_graphics
show_imported_graphic('Bechr29_dual_mountains.png', 30, ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "BeChR29vehdrgMtn", ...
"Effect of GBR-12909 on the reward-mountain data obtained from rat BeChR29");
These shifts can be most clearly discerned in the contour-plot display and bargraph summary:
if show_graphics
show_imported_graphic('BeChR29_quad_GBR.png', 25, ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "BeChR29quad", ...
"Effect of GBR-12909 on the reward-mountain data obtained from rat BeChR29");
This pattern of shifts is seen in the results from 6 of the 7 rats:
if show_graphics
show_imported_graphic('optoGBR_shift_summary.png', 100, ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "GBRshiftSummary", ...
"Effect of GBR-12909 on the reward-mountain data obtained from all 7 rats");
The error bars surrounding the data point for each rat represent the bootstrapped 95% confidence interval (CI) for the parameter in question, whereas the boxes denoted by the heavy solid lines represent the inter-quartile range (IQR) of the median parameter estimate for the group of seven rats. The 
and 
values have been corrected for differential frequency-following fidelity in the vehicle and drug conditions. (See below and main body of the manuscript.)
and values have been corrected for differential frequency-following fidelity in the vehicle and drug conditions. (See below and main body of the manuscript.)Modeling the results
We now undertake modeling to determine how the optical pulse train can be translated into the observed behavioral responses in a manner consistent with the observed effect of dopamine-transporter blockade on the reward mountain.
The first step is to relate the induced frequency of firing in the dopamine neurons to the optical pulse frequency. For this purpose, we adopt a filter function of the form defined by
disp(string({strcat({'Function '},num2str(fun_tab.Number(fun_tab.Name=='FilterFun')))}));
, which provides an accurate description of how medial forebrain bundle fibers subserving electrical intracranial self-stimulation respond to the electrical pulse frequency. It is not known how well the functional form describes frequency-following in optically activated, ChR2-expressing, midbrain dopamine neurons. As explained in the main body of the manuscript, parameter values have been chosen on the basis of published electrophysiological, voltammetric, and behavioral data.
The maximum firing frequency in the dopamine neurons is far lower than in the directly stimulated neurons responsible for the rewarding effect of electrical medial forebrain bundle stimulation (Solomon et al., 2015). To reflect this in the equation for the frequency response, we used updated names and values for the parameters:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "FreqFolDA", "Frequency of firing as a function of pulse frequency");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "FdaBend", "parameter determining sharpness of bend in FreqFol function");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Fopt", "optical pulse frequency");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "FdaRate", "optically induced firing rate of the dopamine neurons");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "FdaRO", "the pulse frequency in the center of the roll-off region");
In later stages of the modeling, we use the inverse function, FilterFunBS, to return an optical pulse frequency given the dopamine firing frequency it produces. (see Functional building blocks for the simulations).
The following graph show the function employed to model frequency-following fidelity in the ChR2-expressing dopamine neurons in response to the optical pulse frequency:
FdaBend = 20;
FdaRO = 50;
logFopt = 0:0.025:2.4;
Fopt = 10.^logFopt;
FRda = FilterFun(Fopt, FdaBend, FdaRO); % Compute firing rate using the frequency roll-off function
FF_graphDA = plot_freqFoll(Fopt,FRda,'DA',1,250,1,100);
if show_graphics
FF_graphDA.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "FreqFolDA", ...
"Induced firing frequency in dopamine neurons as a function of pulse frequency");
The spike counter
Performance for optical stimulation of midbrain dopamine neurons depends both on the induced firing frequency and the number of dopamine neurons excited by the optical input (Ilango et al., 2014). The simplest assumption consistent with this finding is that the behavioral effects depend on the aggregate rate of induced firing, as is the case of the behavioral effects produced by rewarding electrical stimulation of the medial forebrain bundle. The aggregate rate of firing induced by the optical stimulation in the population of dopamine neurons is the product of the number of optically activated neurons and the induced firing frequency:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "AggDAfr", "Aggregate firing rate in the dopamine neurons");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "FdaAgg", "aggregate firing rate of the optically excited dopamine neurons");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Nda", "number of optically excited dopamine neurons");
clear a Fopt FRda logFopt;
clear -regexp ^Fda ^FF_graph;
toc
The effect on the reward mountain produced by modulation of dopaminergic neurotransmission
tic;
close all;
The influence of the induced firing of the dopamine neurons can be modulated by drugs, such as transporter blockers, that alter synaptic transmission. To capture such effects, we add a scalar at the output of the dopamine spike count:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "DAdrive", "scaled post-synaptic influence of the induced dopamine firing");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Kda", "constant that scales the postsynaptic impact of the dopamine firing");
In the following schema, the aggregate rate of firing is represented by the Π symbol and the dopamine-drive scalar by a triangle:
if show_graphics
show_imported_graphic('counter_opto.png',60,ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "DAdrive", ...
"Influence of optically stimulated dopamine release on downstream targets");
We refer to the scaled output of the activated population of midbrain dopamine neurons as "dopamine drive."
The reward-growth function for oICSS
In six of the seven rats, dopamine-transporter blockade shifted the mountain leftward along the pulse-frequency axis. As we demonstrated above, such shifts require that the function that translates dopamine drive into reward intensity (i.e., the reward-growth function) must have a position parameter that is independent of the parameter that sets the maximum reward intensity attainable. The logistic function described by Gallistel's group in the case of eICSS has this property, and we have therefore added such a function at the output of the counter.
if show_graphics
show_imported_graphic('DA_logisitc_simplified.png',50,ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "oICSS_simplified", ...
"simplied oICSS model with logistic reward growth");
In the case of eICSS, the form of the two-dimensional reward-growth function employed in the simulations is based on the empirical measurements obtained by Gallistel's group. No such measurements have yet been obtained for oICSS. Future experiments will be required to determine whether the logistic form or another form with independent input- and output-scaling parameters provides the best description of reward growth in the case of oICSS.
In the case of eICSS, the reward-growth function has been generalized to three dimensions {reward intensity, pulse frequency, train duration} on the basis of additional empirical data (Sonnenschein, Conover, & Shizgal, 2003). Both the pulse frequency and train duration have been manipulated in oICSS experiments (Ilango et al., 2014), but not in a way that makes it possible to describe the joint dependence of reward intensity on these two variables. Future experiments will be required to determine this. For present purposes, we will again adopt the function derived from eICSS data for use in simulating oICSS performance.
We define the position parameter of the reward-growth function for oICSS of midbrain dopamine neurons by modifying
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='fH')))}));
as follows:
[eqn_num, eqn_tab] = ...
add_eqn(eqn_num, eqn_tab, "FHda", "strength-duration function for optical trains delivered to dopamine neurons");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Cda", "chronaxie of the strength-duration function for trains delivered to dopamine neurons");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "Dopt", "duration of an optical pulse train");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "RhoDA", "rheobase of the strength-duration function for trains delivered to dopamine neurons");
The function FpulseHMfun
disp(string({strcat({'Function '},num2str(fun_tab.Number(fun_tab.Name=='FpulseHMfun')))}));
(See Functions composing the reward-mountain model below), accepts an argument list of variable length. If is appended to the end of the argument list, then it will enter into the calculation of 
as specified in the above equation. If is omitted, it assumes an implicit value of one.
is appended to the end of the argument list, then it will enter into the calculation of as specified in the above equation. If is omitted, it assumes an implicit value of one.The following code produces estimates of 
for two values of . The lower value (1) represents the effectiveness of dopaminergic neurotransmission in the vehicle condition, whereas the higher value (2) represents the effectiveness of dopaminergic neurotransmission under the influence of dopamine-transporter blockade.
for two values of . The lower value (1) represents the effectiveness of dopaminergic neurotransmission in the vehicle condition, whereas the higher value (2) represents the effectiveness of dopaminergic neurotransmission under the influence of dopamine-transporter blockade.C = 0.473; % median from Sonnenschein et al., 2003
D = 1.0; % train duration for oICSS study
FdaBend = 20;
FdaRO = 50;
Kda = [1,10^.2];
NnDA = 100;
RhoPiDA = 1572; % See the calculation in the Series-circuit section below for D = 1 s
FPmax = 1000;
FdaHMkDA = FpulseHMfun(C, D, FdaBend, FPmax, FdaRO, NnDA, RhoPiDA, Kda)
%FdaHMkDA is a two-element vector
Note that although Kda(1) is half the value of Kda(2), FdaHMkDA(2) is less than half the value of FdaHMkDA(1). The reason for this is that the frequency-response function for the dopamine neurons has already started to roll off by FdaHMkDA(1) (some dopamine neurons fail to fire once per pulse at this pulse frequency). This will be reflected in a greater discrepancy between RnormMax(1) and one than between RnormMax(2) and one.
To compute time allocation, we also require the value of the parameter that positions the mountain along the price axis, . In the initial simulation, we assume that the value of this parameter is unaffected by the blockade of the dopamine transporter.
. In the initial simulation, we assume that the value of this parameter is unaffected by the blockade of the dopamine transporter.dotPhiObj = 1;
Kaa = 1;
Kec = 1;
Krg = 1;
dotRaa = 0.1;
pObj = 1;
PsubBend = 0.5;
PsubMin = 1.82;
FPmax = 1000; % This value is sufficiently above FdaRO to maximize the firing rate
g = 5;
RnormMax = fRbsrNorm(FPmax, FdaBend, FdaHMkDA, FdaRO, g)
% RnormMax is a two-element vector
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RnormMax)
% PobjE is a two-element vector
We can now estimate time allocation and simulate the reward mountains for the vehicle and drug conditions using the full logistic reward-growth model shown in the figure below:
if show_graphics
show_imported_graphic('oICSS_logistic+3Pe_noMasks_v4.png',25,ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "oICSS_logistic", ...
"oICSS model with logistic reward growth");
We will now simulate the output of this model, disregarding, for the time being, any drug-induced changes in dopamine tone .
a = 3;
Fopt = logspace(0,3,121)'; % column variable
Pobj = logspace(0,3,121); % row variable
Tda1 = TAfun(a, Fopt, FdaBend, FdaHMkDA(1), FdaRO, g, Pobj, PobjE(1), ...
PsubBend, PsubMin, RnormMax(1));
Tda2 = TAfun(a, Fopt, FdaBend, FdaHMkDA(2), FdaRO, g, Pobj, PobjE(2), ...
PsubBend, PsubMin, RnormMax(2));
Before plotting the bar graphs, we need to correct the estimates of the location parameters for the differential frequency-following fidelity in the drug and vehicle conditions. The mountain is shifted downwards along the pulse-frequency axis in the drug condition and thus, frequency-following fidelity is better than in the vehicle condition.
As we explain in the main body of the manuscript,
and
FdaHMkDAstar = FilterFun(FdaHMkDA,FdaBend,FdaRO);
PsubE = PsubFun(PobjE,PsubBend,PsubMin);
PsubEstar = PsubE ./ RnormMax;
PobjEstar = PsubBsFun(PsubEstar,PsubBend,PsubMin);
title_str1 = strcat({'Kda = '}, sprintf('%2.2f', Kda(1)));
xmin = 0;
xmax = 2.3;
ymin = 0.3;
ymax = 2.3;
MTNkDA1 = plot_MTN(Fopt, Pobj, Tda1, 'off', 'MTNkDA1', title_str1, ...
graphs2files, FigDir, xmin, xmax, ymin, ymax);
title_str2 = strcat({'Kda = '}, sprintf('%2.2f', Kda(2)));
MTNkDA2 = plot_MTN(Fopt, Pobj, Tda2, 'off', 'MTNkDA2', title_str2, ...
graphs2files, FigDir, xmin, xmax, ymin, ymax);
MTNkDA1vskDA2 = dual_subplot(MTNkDA1, MTNkDA2, 'MTNkDA1vskDA2',...
graphs2files,FigDir);
if show_graphics
MTNkDA1vskDA2.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "kDA1vskDA2_mtns", ...
"Effect of dopamine-transporter blockade on the reward mountain");
ContkDA1 = plot_contour(Fopt, Pobj, Tda1, PobjE(1), FdaHMkDA(1), 'off', 'ContkDA1', title_str1, ...
strcat({'Kda = '}, num2str(Kda(1))), graphs2files, FigDir,...
xmin, xmax, ymin, ymax);
ContkDA2 = plot_contour(Fopt, Pobj, Tda2, PobjE(2), FdaHMkDA(2), 'off', 'ContkDA2', title_str2, ...
strcat({'Kda = '}, num2str(Kda(2))), graphs2files, FigDir,...
xmin, xmax, ymin, ymax);
bg_kDA1vskDA2 = plot_bgStar(FdaHMkDA(1), FdaHMkDA(2), FdaHMkDAstar(1), FdaHMkDAstar(2),...
PobjE(1), PobjE(2), PobjEstar(1), PobjEstar(2),...
'off', 'kDA1vskDA2_bg', graphs2files, FigDir);
bg_root = 'bg_kDA1vskDA2';
quad_kDA1vskDA2 = quad_subplot(ContkDA1, ContkDA2, bg_kDA1vskDA2, 'quad_kDA1vskDA2', bg_root, ...
graphs2files, FigDir);
if show_graphics
quad_kDA1vskDA2.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "kDA1vskDA2_quad", ...
"Effect of dopamine-transporter blockade on the reward mountain");
In the above simulation, the drug-induced boost in the effectiveness of dopaminergic neurotransmission has rescaled the input to the reward-growth function. As a result, the mountain shifts leftwards along the pulse-frequency axis. This captures one of the observed effects of GBR-12909.
The small shift along the price axis is due to the fact that is closer to one in the drug condition than in the vehicle condition. This is due to the fact that is lower in the drug condition and thus further from the pulse frequency at which frequency following begins to roll off. The corrected value of 
is designated by the dot-dash cyan line on the red bar.
is closer to one in the drug condition than in the vehicle condition. This is due to the fact that is lower in the drug condition and thus further from the pulse frequency at which frequency following begins to roll off. The corrected value of is designated by the dot-dash cyan line on the red bar.The difference in frequency-following fidelity in the drug and vehicle conditions also produces a small shift in 
. This shift is removed by the correction (dot-dash cyan line on the blue bar).
. This shift is removed by the correction (dot-dash cyan line on the blue bar).The empricially observed rightward shifts along the pulse axis survive the correction for differential frequency-following fidelity. To capture these, a scalar increase must be produced at or beyond the output of the reward-growth function. In the example shown below, we do this by having the drug-induced increase in dopamine tone reduce the subjective effort cost, as proposed by Salomone and colleagues.
clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
Dopaminergic modulation of subjective effort cost
tic;
close all;
In
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='oICSS_logistic')))}));
, blockade of the dopamine transporter changes both phasic and tonic dopamine signaling. Given that the lack of correlation between the magnitudes of the resulting shifts along the pulse-frequency and price axes, separate scalars are included for phasic and tonic signaling. The schema indicates three different ways in which increases in tonic signaling could shift the mountain along the price axis, as designated by the black dashed line. The one simulated here is a reduction in the rate of subjective exertion entailed in holding down the lever.
% Set Fhm to the values obtained from ratBechr29 in the vehicle and drug conditions
Kda = [1,1.4458]; % Kda value for Bechr29
FdaHMkDA = [31.34679397, 20.77560124]; % values for Bechr29
FdaHMkDAstar = [26.29253041, 18.18388501]; % values for Bechr29
PsubBend = 0.5;
PsubMin = 1.8197;
FdaBend = 20;
FdaRO = 50;
FPmax = 200; % This is the value that was employed in the correction of the location parameters
g = 3.160774009; % value for Bechr29; value common for Veh and Drg conditions
RnormMax = [0.893580092, 0.964209143]; % means of resampled values for Bechr29
% PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RnormMax, KeffMod)
PobjE = [6.295449067, 10.04096421]; % values for Bechr29
PobjEstar = [7.375024269, 10.47998071]; % values for Bechr29
% We can't simply perform the corrections of the location parameters on the means.
% The original corrections were performed on the 250 resampled log values.
% The statistics are based on these vectors.
% The log of the mean is not equal to the mean of the logs.
% Similarly, we cannot simply use FilterFun to correct the mean Fhm values and use these to estimate RnormMax
% logFfiringHM (i.e., logFhmStar) and logPobjEstar were estimated from each of the 250 resampled values of
% logFpulseHM and logPobjE. We encounter the same issue.
% The values used above are based on the 250 resampled log values, which were loaded directly from the
% saved workspace that was used in the fits.
a = 3;
Fopt = logspace(0,3,121)'; % column variable
Pobj = logspace(0,3,121); % row variable
Tda1 = TAfun(a, Fopt, FdaBend, FdaHMkDA(1), FdaRO, g, Pobj, PobjE(1), ...
PsubBend, PsubMin, RnormMax(1));
Tda2 = TAfun(a, Fopt, FdaBend, FdaHMkDA(2), FdaRO, g, Pobj, PobjE(2), ...
PsubBend, PsubMin, RnormMax(2));
xmin = 0;
xmax = 2.3;
ymin = 0.3;
ymax = 2.3;
title_str1 = strcat({'Kda = '}, sprintf('%2.2f', Kda(1)));
MTNkDAeffMod1 = plot_MTN(Fopt, Pobj, Tda1, 'off', 'MTNoICSSkDAeffMod1', title_str1, ...
graphs2files, FigDir, xmin, xmax, ymin, ymax);
title_str2 = strcat({'Kda = '}, sprintf('%2.2f', Kda(2)));
MTNkDAeffMod2 = plot_MTN(Fopt, Pobj, Tda2, 'off', 'MTNoICSSkDAeffMod2', title_str2, ...
graphs2files, FigDir, xmin, xmax, ymin, ymax);
MTNoICSSkDA1vskDA2 = dual_subplot(MTNkDAeffMod1, MTNkDAeffMod2, 'MTNoICSSkDA1vskDA2',...
graphs2files,FigDir);
if show_graphics
MTNoICSSkDA1vskDA2.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "kDA1vskDAeffMod2_mtns", ...
"Effect of dopamine-transporter blockade on the reward mountain");
ContkDAeffMod1 = plot_contour(Fopt, Pobj, Tda1, PobjE(1), FdaHMkDA(1), 'off', 'ContkDAeffMod1', title_str1', ...
strcat({'Kda = '}, num2str(Kda(1))), graphs2files, FigDir,...
xmin, xmax, ymin, ymax);
ContkDAeffMod2 = plot_contour(Fopt, Pobj, Tda2, PobjE(2), FdaHMkDA(2), 'off', 'ContkDAeffMod2', title_str2, ...
strcat({'Kda = '}, num2str(Kda(2))), graphs2files, FigDir,...
xmin, xmax, ymin, ymax);
ymin = -0.3581; % from the bargraph for Bechr29
ymax = 0.3652; % from the bargraph for Bechr29
bg_kDA1vskDA2effMod = plot_bgStar(FdaHMkDA(1), FdaHMkDA(2), FdaHMkDAstar(1), FdaHMkDAstar(2),...
PobjE(1), PobjE(2), PobjEstar(1), PobjEstar(2),...
'off', 'kDA1vskDA2_bg', graphs2files, FigDir,...
ymin, ymax);
bg_root = 'bg_kDA1vskDAeffMod2';
quad_kDA1vskDA2effMod = quad_subplot(ContkDAeffMod1, ContkDAeffMod2, bg_kDA1vskDA2effMod,...
'quad_kDA1vskDAeffMod2', bg_root, graphs2files, FigDir);
if show_graphics
quad_kDA1vskDA2effMod.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "kDA1vskDA2effMod_quad", ...
"Effect of dopamine-transporter blockade on the reward mountain");
The simulated results are now in qualitative accord with the empirical findings: the mountain shifts leftwards along the pulse-frequency axis and rightwards along the price axis in response to blockade of the dopamine transporter.
Here are the shifts observed in the mountains obtained from rat BeChR29 and the shifts simulated above using the and 
values derived from that rat's data:
and values derived from that rat's data:BeChR29_bg = make_fig_from_png(fullfile(ImpFigDir,'BeChR29_bg.png'), 30);
BeChR29bgVSsimBG = dual_subplot(BeChR29_bg, bg_kDA1vskDA2effMod, 'BeChR29bgVSsimBG',...
graphs2files,FigDir);
% rescale the right panel to 60% width, 80% height & center in right panel
BeChR29bgVSsimBGrs = adjust_right_panel(BeChR29bgVSsimBG, 0.6, 0.825);
if show_graphics
BeChR29bgVSsimBGrs.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "BeChR29bgVSsimBGrs", ...
"Simulated and observed effects of dopamine-transporter blockade on the reward mountain");
The correspondence is forced by inputting values of 
, 
, and g obtained from the fit of the reward-mountain model to the empirical data. All that this comparison does is to verify the functions for computing time allocation and graphing the results.
, , and g obtained from the fit of the reward-mountain model to the empirical data. All that this comparison does is to verify the functions for computing time allocation and graphing the results.Full simulation would yield the same results if we adjusted 
so as to generate the observed values.
so as to generate the observed values.clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
The series-circuit model of eICSS and oICSS
tic;
close all;
In the previous section, we demonstrate how the effects of dopamine-transporter blockade on oICSS can be explained within the context of the mountain model by a combination of changes in phasic and tonic dopamine signaling. In this section, we pose an additional challenge: can the explanation proposed for the oICSS data also account for the known effects of dopamine-transporter blockade on eICSS? Below, we show that the proposed explanation fails to meet this challenge, we discuss the implications for the "series-circuit" model of eICSS, and we develop a new model that can account for the effects of dopamine-transporter blockade on both oICSS and eICSS.
The series-circuit model treats oICSS and eICSS as behavioral manifestations of the effects produced by injecting signals at two different neural stages of the same pathway. On this view, eICSS arises from electrically induced activation of highly excitable, non-dopaminergic neurons that project directly or indirectly to midbrain dopamine neurons. In other words, the directly activated neurons that give rise to eICSS are in series with the dopamine neurons that render the electrical stimulation rewarding. The midbrain dopamine neurons are excited trans-synaptically in the case of eICSS and directly in the case of oICSS. The consequences of their activation are the same: the subject seeks to re-initiate the stimulation and will pay large effort and opportunity costs when the strength of the stimulation suffices to produce a large increment in the aggregate firing rate of the dopamine neurons.
Here, we show only the portion of the model that generates the reward-intensity signal.
if show_graphics
show_imported_graphic('series-circuit_DA_double_logisitc_RG_v2.png',15,ImpFigDir);
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "SeriesCircuit_eICSS_oICSS", ...
"Series-circuit model of eICSS and oICSS");
As the figure above shows, two different reward-growth functions are required, one upstream of the dopamine neurons and another downstream. The upstream reward-growth function is necessary in order to account for the data summarized in
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='PimGBRshifts')))}));
That figure shows that modulation of dopaminergic neurotransmission alters eICSS performance by shifting the reward mountain along the price axis and not along the pulse-frequency axis. This implies that the drug-induced change in dopamine signaling acts at or beyond the output of the reward-growth function for eICSS. For this to occur in the series-circuit model, the reward-growth function for eICSS must lie upstream (to the left) of the dopamine neurons.
A second reward-growth function must lie downstream (to the right) of the dopamine neurons. This is required because GBR-12909 shifted the reward mountain for oICSS along the pulse-frequency axis, as shown in
disp(string({strcat({'Figures '},num2str(fig_tab.Number(fig_tab.Name=='BeChR29vehdrgMtn')))}));
and
disp(string({strcat({''},num2str(fig_tab.Number(fig_tab.Name=='GBRshiftSummary')))}));
We explain above (The significance of orthogonal shifts) that the mountain will move along the pulse-frequency axis only if the input to the reward-growth function has been rescaled. Such rescaling occurs when the magnitude of dopamine transients is boosted by GBR-12909 and the reward-growth function for oICSS is positioned downstream of the dopamine neurons.
Series-circuit model: changes in the reward mountain for eICSS in response to dopamine-transporter blockade
To generate a reward mountain for eICSS from the series-circuit model, we proceed in several stages. As was done above (The reward-growth function for oICSS), the optical pulse frequency required to produce a half-maximal reward intensity is calculated directly from
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='fH')))}));
We next need to compute the electrical pulse frequency (applied to a medial forebrain bundle electrode) that produces excitation in the dopamine neurons equivalent to that produced by a given train of optical pulses delivered to the midbrain dopamine neurons. According the series-circuit model, all inputs that produce the same peak output from the dopamine neurons will produce the same rewarding effect. This will be true regardless of whether the dopamine neurons are excited directly by optical activation or indirectly by trans-synaptic input from medial forebrain bundle neurons activated by electrical stimulation. According to the model, an observer placed downstream from the dopamine neurons and supplied only with information about the aggregate peak output of these neurons cannot know whether optical or electrical stimulation was responsible for a given phasic increase in dopamine release. It is the peak magnitude of this phasic increase in aggregate firing that determines the intensity of the rewarding effect.
(See "The spike counter" above.)
To obtain the electrical pulse frequency required to produce a reward of half-maximal intensity, we need to back-solve the equations describing the stages of the model that intervene between the electrode and the dopamine neurons. The back-solutions return the electrical pulse frequency that delivers an input to dopamine neurons equivalent to the optical pulse frequency that produces a half-maximal reward intensity.
The reward-growth functions (S-shaped curves in rectangular boxes) in the flow diagrams are normalized: Their output varies from zero to one and is then scaled by the variable in the triangle to their right. In the case of the upstream reward-growth function in
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='SeriesCircuit_eICSS_oICSS')))}));
(the reward-growth function to the left of the dopamine neurons), that scaling variable is 
. The value of this variable determines the maximum input that the electrode can deliver to the dopamine neurons, scaled in terms of the equivalent optical pulse frequency. In the initial simulation below, is set to 63 pulses per second, and given the parameters of the frequency-following function of the dopamine neurons {
}, this will drive the dopamine neurons to a near-maximal level (~43 spikes ). We demonstrate below that the qualitative effect of dopamine transporter blockade is little affected by the value of .
. The value of this variable determines the maximum input that the electrode can deliver to the dopamine neurons, scaled in terms of the equivalent optical pulse frequency. In the initial simulation below, is set to 63 pulses per second, and given the parameters of the frequency-following function of the dopamine neurons {}, this will drive the dopamine neurons to a near-maximal level (~43 spikes ). We demonstrate below that the qualitative effect of dopamine transporter blockade is little affected by the value of .The function that computes the required electrical pulse frequency in the series-circuit model, FscHMbs, comprises three back-solutions. (See Functions composing the reward-mountain model below.) First, we invert the scaling of the output of the upstream reward-growth function: we divide the optical pulse frequency required to produce a reward of half-maximal intensity by . This gives us the output of the upstream normalized reward-growth function (a value between zero and one), the fraction of required to produce a half-maximal reward intensity at the output of the dopamine neurons (the FdaHMkDA value). We then backsolve the upstream reward-growth function (by means of the LogistNormBsFun function) to obtain the average firing rate of the medial forebrain bundle neurons required to produce the reward intensity in question, given the position parameter of the logistic reward-growth function and the value of its exponent. The position parameter of the logistic is calculated using
. This gives us the output of the upstream normalized reward-growth function (a value between zero and one), the fraction of required to produce a half-maximal reward intensity at the output of the dopamine neurons (the FdaHMkDA value). We then backsolve the upstream reward-growth function (by means of the LogistNormBsFun function) to obtain the average firing rate of the medial forebrain bundle neurons required to produce the reward intensity in question, given the position parameter of the logistic reward-growth function and the value of its exponent. The position parameter of the logistic is calculated using disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='fH')))}));
(the FFmfbHM function). Finally, we use the FilterFunBS function to obtain the electrical pulse frequency that produces this average firing rate.
The input to FscHMbs is the optical pulse frequency required to produce a reward of half-maximal intensity (
) by direct activation of the dopamine neurons, and the output is the equivalent electrical pulse frequency (
).
) by direct activation of the dopamine neurons, and the output is the equivalent electrical pulse frequency ().[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "FscHMbs", ...
"FhmDA, FmfbBend, FhmMFB, FmfbRO, gMFB, KrgEq",...
"Back solution of the full reward-growth function for the series-circuit (sc) model");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "FscHM", "position parameter of RG function in the series-circuit model");
We now compute the value of the FscHMbs function for the vehicle and drug conditions.
%% Calculate FscHM for the vehicle and drug conditions
% Calculate the position parameter of the upstream reward-growth function
C = 0.473; % median from Sonnenschein et al., 2003
D = 0.5; % typical train duration
FdaBend = 20;
logFmfbBend = 1.3222; % from Solomon et al., 2015
FmfbBend = 10^logFmfbBend
FmfbBend = 20.9991
FdaRO = 50;
gMFB = 1.58; % reduce the value due to the embedding within the DA RG function
logFmfbRO = 2.5587; % from Solomon et al., 2015
FmfbRO = 10^logFmfbRO
FmfbRO = 361.9929
NnMFB = 126;
RhoPiMFB = 5000; % RhoPi/N is within ~20% of 50, which is roughly consistent with Sonnenchein et al., 2003
FPmax = 1000;
FmfbHM = FpulseHMfun(C, D, FmfbBend, FPmax, FmfbRO, NnMFB, RhoPiMFB)
FmfbHM = 77.2222
In the vehicle condition of the current experiment, the median value for oICSS was 27.1. The train duration was 1 s, whereas it was 0.5 s in the corresponding eICSS study. The form and parameters of the temporal-integration function for oICSS of midbrain dopamine neurons are unknown. Faut de mieux, we will use the functional form and parameters obtained from eICSS to compute the Fhm value for a 0.5 s train that corresponds to the value obtained in the current oICSS study for 1 s trains.
value for oICSS was 27.1. The train duration was 1 s, whereas it was 0.5 s in the corresponding eICSS study. The form and parameters of the temporal-integration function for oICSS of midbrain dopamine neurons are unknown. Faut de mieux, we will use the functional form and parameters obtained from eICSS to compute the Fhm value for a 0.5 s train that corresponds to the value obtained in the current oICSS study for 1 s trains. In the series-circuit model, the output of the MFB neurons must pass through the midbrain dopamine neurons. If so, the temporal-integration characteristics obtained in eICSS experiments reflect those of both the directly stimulated and dopamine stages of the circuit. According to this model, integration in the dopamine stage cannot be faster than estimated in the eICSS study.
It follows from
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='fH')))}));
that the Fhm value for a 0.5 s train that corresponds to an Fhm value for a 1 s train is given by
[eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, ...
"FfiringD1D2", "calculates equivalently effective pulse frequencies at two different train durations");
Given the frequency-following function, the median value for the vehicle condition (27.1) corresponds to a firing frequency of 23.1.
value for the vehicle condition (27.1) corresponds to a firing frequency of 23.1. FpulseHMmedVeh = 27.094;
FfiringHMmedVeh = FilterFun(FpulseHMmedVeh, FdaBend, FdaRO)
FfiringHMmedVeh = 23.1475
Given these values
C = 0.473;
D1 = 1;
D2 = 0.5;
FfiringHmD1 = FfiringHMmedVeh
FfiringHmD1 = 23.1475
FfiringHmD2 = FfiringHmD1 .* (1+(C./D2)) ./ (1+(C./D1))
FfiringHmD2 = 30.5805
FPmax = 1000;
FpulseHmD2 = FilterFunBS(FdaBend, FdaRO, FfiringHmD2,FPmax)
FpulseHmD2 = 37.6179
The MFB will have to deliver excitation equivalent to 
to deliver a half-maximal reward intensity from the dopamine neurons at a train duration of 0.5 s. We will determine the values of and 
so as to obtain the above 
value in the vehicle condition.
to deliver a half-maximal reward intensity from the dopamine neurons at a train duration of 0.5 s. We will determine the values of and so as to obtain the above value in the vehicle condition. We first set the number of activated dopamine neurons arbitrarily to 100.
NnDA = 100;
Then, we use
disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='FfiringD1D2')))}));
to find the required value of , keeping in mind that is the aggregate firing rate required to produce a half-maximal reward intensity with a train of infinite duration. Setting D1 to infinity makes the denominator of
, keeping in mind that is the aggregate firing rate required to produce a half-maximal reward intensity with a train of infinite duration. Setting D1 to infinity makes the denominator of disp(string({strcat({'Equation '},num2str(eqn_tab.Number(eqn_tab.Name=='FfiringD1D2')))}));
equal to one, so the equation is simplified to:
RhoPiDA = NnDA * (FfiringHmD2/(1+(C/D2)))
RhoPiDA = 1.5715e+03
We are now positioned to compute the MFB pulse frequency that will drive the dopamine neurons to produce a reward of half-maximal intensity.
Kda = [1,10^0.15]; % Effect of dopamine-transporter blockade on the frequency of firing required
% to produce a reward of half-maximal intensity. Median Kda in current oICSS study: 10^0.1446
FPmax = 1000; % This pulse frequency well above FdaRO and can thus serve to estimate their firing rate
FdaHMkDA = FpulseHMfun(C, D2, FdaBend, FPmax, FdaRO, NnDA, RhoPiDA, Kda) % result is a two-element vector
Now, we work backwards through the MFB input to find the electrical pulse frequency that will drive the dopamine neurons to produce a reward of half-maximal intensity.
KrgEq = 63;
FFmax = FilterFun(FPmax, FmfbBend, FmfbRO);
if FFmax > LogistNormBsFun(gMFB,FmfbHM,max(FdaHMkDA)/KrgEq)
for j = 1:length(FdaHMkDA)
FscHM(j) = FscHMbs(FdaHMkDA(j), FmfbBend, FmfbHM, FmfbRO, FPmax, gMFB, KrgEq) % result is a two-element vector
end
else
display("FscHM could not be calculated because KrgEq is too low.");
display("Set KrgEq such that FFmax > FilterFun(FPmax, Fbend, Fro).");
return
end
FscHM = 99.0574
% Calculate PobjE for the vehicle and drug conditions
dotPhiObj = 1;
Kaa = 1;
Kec = 1;
KeffMod = Kda; % We tie the subjective effort cost to the effect of the drug
Kda
% Median shift in Pe in the current oICSS study: 10^0.1525
Krg = 1;
dotRaa = 0.1;
pObj = 1;
PsubBend = 0.5;
PsubMin = 1.82;
Felec = logspace(0,3,121)'; % column variable
To solve the time-allocation equation for the series-circuit model, we must first compute the drive on the dopamine neurons that is produced by each value of the electrical pulse frequency (
). We express this drive in terms of the optical pulse frequency that produces equivalent firing in the dopamine neurons, 
. The FoptEquivFun function first translates the electrical pulse frequency into the induced rate of firing in the medial forebrain bundle neurons, then translates this firing rate into a normalized reward intensity, and then converts the normalized reward intensity into the equivalent optical pulse frequency. Once this has been done, time allocation can be computed using the same equation and parameter values that were used to generate the reward mountain for optical stimulation. The FoptEquivFun function is the inverse of the FscHMbs function. (See Functions composing the reward-mountain model below.)
). We express this drive in terms of the optical pulse frequency that produces equivalent firing in the dopamine neurons, . The FoptEquivFun function first translates the electrical pulse frequency into the induced rate of firing in the medial forebrain bundle neurons, then translates this firing rate into a normalized reward intensity, and then converts the normalized reward intensity into the equivalent optical pulse frequency. Once this has been done, time allocation can be computed using the same equation and parameter values that were used to generate the reward mountain for optical stimulation. The FoptEquivFun function is the inverse of the FscHMbs function. (See Functions composing the reward-mountain model below.) [fun_num, fun_tab] = add_fun(fun_num, fun_tab, "FoptEquivFun", ...
"Felec, FmfbBend, FmfbHM, FmfbRO, gMFB, KrgE",...
"Function to compute the optical pulse frequency that produces the same DA output as Felec");
[sym_num, sym_tab] = add_sym(sym_num, sym_tab, "FoptEquiv", "optical pulse frequency that applies the same drive on the dopamine neurons as Felec");
FoptEquiv = FoptEquivFun(Felec, FmfbBend, FmfbHM, FmfbRO, gMFB, KrgEq);
FPmaxEq = max(FoptEquiv) % This is the highest equivalent optical pulse frequency that the DA neurons will "see"
FPmaxEq = 57.9538
gDA = 5; % Use a value consistent with the oICSS study; median in vehicle condition was 4.61
RnormMax = fRbsrNorm(FPmaxEq, FdaBend, FdaHMkDA, FdaRO, gDA)
% RnormMax is a 2-element vector
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RnormMax, KeffMod)
% PobjE is a two-element vector
Pobj = logspace(0,3,121); % row variable
% Compute time allocation for the vehicle and drug conditions
a = 2; % median in vehicle condition was 1.84
Tsc1 = TAfun(a, FoptEquiv, FdaBend, FdaHMkDA(1), FdaRO, gDA, Pobj, ...
PobjE(1), PsubBend, PsubMin, RnormMax(1));
Tsc2 = TAfun(a, FoptEquiv, FdaBend, FdaHMkDA(2), FdaRO, gDA, Pobj, ...
PobjE(2), PsubBend, PsubMin, RnormMax(2));
title_str1 = strcat({'Kda = '}, sprintf('%2.2f', Kda(1)));
MTNkDA1sc = plot_MTN(Felec, Pobj, Tsc1, 'off', 'MTNkDA1sc', title_str1, ...
graphs2files, FigDir);
title_str2 = strcat({'Kda = '}, sprintf('%2.2f', Kda(2)));
MTNkDA2sc = plot_MTN(Felec, Pobj, Tsc2, 'off', 'MTNkDA2sc', title_str2, ...
graphs2files, FigDir);
MTNkDA1vskDA2sc = dual_subplot(MTNkDA1sc, MTNkDA2sc, 'MTNkDA1vskDA2sc',...
graphs2files,FigDir);
if show_graphics
MTNkDA1vskDA2sc.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "kDA1vskDA2_sc_mtns", ...
"Effect of dopamine-transporter blockade on the reward mountain");
FscHMstar = FilterFun(FdaHMkDA,FdaBend,FdaRO);
% Computed on the basis of frequency following in the dopamine neurons
PsubE = PsubFun(PobjE,PsubBend,PsubMin);
PsubEstar = PsubE ./ RnormMax;
PobjEstar = PsubBsFun(PsubEstar,PsubBend,PsubMin);
ContkDA1sc = plot_contour(Felec, Pobj, Tsc1, PobjE(1), FscHM(1), 'off', 'ContkDA1sc', title_str1, ...
strcat({'Kda = '}, num2str(Kda(1))), graphs2files, FigDir);
ContkDA2sc = plot_contour(Felec, Pobj, Tsc2, PobjE(2), FscHM(2), 'off', 'ContkDA2sc', title_str2, ...
strcat({'Kda = '}, num2str(Kda(2))), graphs2files, FigDir);
bg_kDA1vskDA2sc = plot_bgStar(FscHM(1), FscHM(2), FscHMstar(1), FscHMstar(2),...
PobjE(1), PobjE(2), PobjEstar(1), PobjEstar(2),...
'off', 'kDA1vskDA2sc_bg', graphs2files, FigDir);
bg_root = 'bg_kDA1vskDA2sc';
quad_kDA1vskDA2sc = quad_subplot(ContkDA1sc, ContkDA2sc, bg_kDA1vskDA2sc, 'quad_kDA1vskDA2sc', bg_root, ...
graphs2files, FigDir);
if show_graphics
quad_kDA1vskDA2sc.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "kDA1vskDA2_sc_quad", ...
"Effect of dopamine-transporter blockade on the reward mountain");
toc
tic
clear -regexp ^Cont ^dual ^MTN ^quad;
The following two figures compare the shifts in the position of the reward mountain predicted by the series-circuit model and observed empirically when eICSS is challenged with GBR-12909, a dopamine-transporter blocker:
GBR_eICSS_summary = make_fig_from_png(fullfile(ImpFigDir,'GBR_eICSS_summary.png'), 100);
GBReICSSVSsimBG = dual_subplot(GBR_eICSS_summary, bg_kDA1vskDA2sc, 'GBReICSSVSsimBG',...
graphs2files,FigDir);
% rescale the right panel to 60% width, 80% height & center in right panel
GBReICSSVSsimBGrs = adjust_right_panel(GBReICSSVSsimBG, 0.75, 1);
if show_graphics
GBReICSSVSsimBGrs.Visible = 'on';
end
The predicted shift along the price axis is consistent qualitatively with the empirical results of the experiment in which eICSS was challenged with GBR-12909: the mountain shifts rightwards along the price axis. However, the predicted shift along the pulse-frequency axis is inconsistent with the empirical results. Whereas the series-circuit model predicts a leftward shift along the pulse-frequency axis, the mountain was not displaced systematically along the pulse-frequency axis in the empirical study. This discrepancy reveals a serious flaw in the series-circuit model.
The reason for the failure of the series-circuit model becomes clear upon inspection of
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='SeriesCircuit_eICSS_oICSS')))}));
The change in dopaminergic neurotransmission produced by GBR-12909 rescales the input to the downstream reward-growth function (to the right of the dopamine neurons), changing its position parameter. In contrast to what is observed empirically in the eICSS studies, this causes the simulated mountain to shift along the pulse-frequency axis, regardless of whether the dopamine neurons are activated directly by optical stimulation or indirectly by synaptic input from the medial forebrain bundle neurons stimulated by an electrode.
An observer positioned at the input to the downstream reward-growth function cannot know whether optical or electrical stimulation was responsible for the phasic dopamine signal that constitutes the input to this function. Thus, the series-circuit model cannot readily generate differential predictions in response to optical and electrical inputs.
toc
The qualitative predictions don't depend meaningfully on the value of the parameter than scales the medial forebrain bundle drive on the dopamine neurons.
tic;
clear -regexp ^GBR;
if show_graphics
close all;
end
The next set of graphs shows that the predictions don't change qualitatively when the maximum medial forebrain bundle drive on the dopamine neurons is reduced: The mountain continues to shift along both the price and pulse-frequency axes
KrgEqLo = 41; % Maximum MFB drive is reduced
% n.b. KrgEqLo must be sufficiently exceed FdaHMkDA so as to keep FscHMlo below the maximum firing rate
% Otherwise the back-solution function will generate an error message and return.
% Re-calculate the location-parameter values along the pulse-frequency axis for
% the series-circuit model as a whole.
FFmax = FilterFun(FPmax, FmfbBend, FmfbRO);
if FFmax > LogistNormBsFun(gMFB,FmfbHM,max(FdaHMkDA)/KrgEqLo)
for j=1:length(FdaHMkDA)
FscHMlo(j) = FscHMbs(FdaHMkDA(j), FmfbBend, FmfbHM, FmfbRO, FPmax, gMFB, KrgEqLo) % result is a two-element vector
end
else
display("FscHM could not be calculated because KrgEqLo is too low.");
display("Set KrgEqLo such that FFmax > FilterFun(FPmax, Fbend, Fro).");
return
end
FscHMlo = 380.5649
% Recalculate the drive on the DA neurons
FoptEquivLo = FoptEquivFun(Felec, FmfbBend, FmfbHM, FmfbRO, gMFB, KrgEqLo);
FPmaxEq = max(FoptEquivLo) % This is the highest equivalent optical pulse frequency that the DA neurons will "see"
FPmaxEq = 37.7159
RnormMax = fRbsrNorm(FPmaxEq, FdaBend, FdaHMkDA, FdaRO, gDA)
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RnormMax, KeffMod)
% Compute time allocation for the vehicle and drug conditions
Tsc3 = TAfun(a, FoptEquivLo, FdaBend, FdaHMkDA(1), FdaRO, gDA, Pobj, ...
PobjE(1), PsubBend, PsubMin, RnormMax(1));
Tsc4 = TAfun(a, FoptEquivLo, FdaBend, FdaHMkDA(2), FdaRO, gDA, Pobj, ...
PobjE(2), PsubBend, PsubMin, RnormMax(2));
title_str1 = strcat({'Kda = '}, sprintf('%2.2f', Kda(1)));
MTNkDA1scKrgEqLo = plot_MTN(Felec, Pobj, Tsc3, 'off', 'MTNkDA1scKrgEqLo', title_str1, ...
graphs2files, FigDir);
title_str2 = strcat({'Kda = '}, sprintf('%2.2f', Kda(2)));
MTNkDA2scKrgEqLo = plot_MTN(Felec, Pobj, Tsc4, 'off', 'MTNkDA2scKrgEqLo', title_str2, ...
graphs2files, FigDir);
MTNkDA1vskDA2scKrgEqLo = dual_subplot(MTNkDA1scKrgEqLo, MTNkDA2scKrgEqLo, 'MTNkDA1vskDA2scKrgEqLo',...
graphs2files,FigDir);
if show_graphics
MTNkDA1vskDA2scKrgEqLo.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "kDA1vskDA2_scKrgEqLo_mtns", ...
"Effect of dopamine-transporter blockade on the reward mountain");
FscHMstar = FilterFun(FdaHMkDA,FdaBend,FdaRO);
% Computed on the basis of frequency following in the dopamine neurons
PsubE = PsubFun(PobjE,PsubBend,PsubMin);
PsubEstar = PsubE ./ RnormMax;
PobjEstar = PsubBsFun(PsubEstar,PsubBend,PsubMin);
ContkDA1scKrgEqLo = plot_contour(Felec, Pobj, Tsc3, PobjE(1), FscHMlo(1), 'off', 'ContkDA1scKrgEqLo', title_str1, ...
strcat({'Kda = '}, num2str(Kda(1))), graphs2files, FigDir);
ContkDA2scKrgEqLo = plot_contour(Felec, Pobj, Tsc4, PobjE(2), FscHMlo(2), 'off', 'ContkDA2scKrgEqLo', title_str2, ...
strcat({'Kda = '}, num2str(Kda(2))), graphs2files, FigDir);
bg_kDA1vskDA2scKrgEqLo = plot_bgStar(FscHM(1), FscHM(2), FscHMstar(1), FscHMstar(2),...
PobjE(1), PobjE(2), PobjEstar(1), PobjEstar(2),...
'off', 'kDA1vskDA2scKrgEqLo_bg', graphs2files, FigDir);
bg_root = 'bg_kDA1vskDA2scKrgEqLo';
quad_kDA1vskDA2scKrgEqLo = quad_subplot(ContkDA1scKrgEqLo, ContkDA2scKrgEqLo, bg_kDA1vskDA2scKrgEqLo, ...
'quad_kDA1vskDA2scKrgEqLo', bg_root, graphs2files, FigDir);
if show_graphics
quad_kDA1vskDA2scKrgEqLo.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "kDA1vskDA2_scKrgEqLo_quad", ...
"Effect of dopamine-transporter blockade on the reward mountain");
toc
tic;
clear -regexp ^bg ^Cont ^dual ^MTN ^quad;
close all;
Next, we increase the medial forebrain bundle drive so that some of it is wasted. The firing rate of the dopamine neurons now asymptotes at sub-maximal levels of medial forebrain bundle drive. Again, the predictions of the series-circuit model do not change qualitatively: The mountain continues to shift along both the price and pulse-frequency axes.
KrgEqHi = 126; % Maximum MFB drive exceeds the DA frequency-following range
% Re-calculate the location-parameter values along the pulse-frequency axis for
% the series-circuit model as a whole.
FFmax = FilterFun(FPmax, FmfbBend, FmfbRO);
if FFmax > LogistNormBsFun(gMFB,FmfbHM,max(FdaHMkDA)/KrgEqHi)
for j=1:length(FdaHMkDA)
FscHMhi(j) = FscHMbs(FdaHMkDA(j), FmfbBend, FmfbHM, FmfbRO, FPmax, gMFB, KrgEqHi) % result is a two-element vector
end
else
display("FscHM could not be calculated because KrgEqHi is too low.");
display("Set KrgEqHi such that FFmax > FilterFun(FPmax, Fbend, Fro).");
return
end
FscHMhi = 44.9730
% Recalculate the drive on the DA neurons
FoptEquivHi = FoptEquivFun(Felec, FmfbBend, FmfbHM, FmfbRO, gMFB, KrgEqHi);
FPmaxEq = max(FoptEquivHi) % This is the highest equivalent optical pulse frequency that the DA neurons will "see"
FPmaxEq = 115.9075
RnormMax = fRbsrNorm(FPmaxEq, FdaBend, FdaHMkDA, FdaRO, gDA)
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RnormMax, KeffMod)
% Compute time allocation for the vehicle and drug conditions
Tsc5 = TAfun(a, FoptEquivHi, FdaBend, FdaHMkDA(1), FdaRO, gDA, Pobj, ...
PobjE(1), PsubBend, PsubMin, RnormMax(1));
Tsc6 = TAfun(a, FoptEquivHi, FdaBend, FdaHMkDA(2), FdaRO, gDA, Pobj, ...
PobjE(2), PsubBend, PsubMin, RnormMax(2));
title_str1 = strcat({'Kda = '}, sprintf('%2.2f', Kda(1)));
MTNkDA1scKrgEqHi = plot_MTN(Felec, Pobj, Tsc5, 'off', 'MTNkDA1scKrgEqHi', title_str1, ...
graphs2files, FigDir);
title_str2 = strcat({'Kda = '}, sprintf('%2.2f', Kda(2)));
MTNkDA2scKrgEqHi = plot_MTN(Felec, Pobj, Tsc6, 'off', 'MTNkDA2scKrgEqHi', title_str2, ...
graphs2files, FigDir);
MTNkDA1vskDA2scKrgEqHi = dual_subplot(MTNkDA1scKrgEqHi, MTNkDA2scKrgEqHi, 'MTNkDA1vskDA2scKrgEqHi',...
graphs2files,FigDir);
if show_graphics
MTNkDA1vskDA2scKrgEqHi.Visible = 'on';
end
FscHMstar = FilterFun(FdaHMkDA,FdaBend,FdaRO);
% Computed on the basis of frequency following in the dopamine neurons
PsubE = PsubFun(PobjE,PsubBend,PsubMin);
PsubEstar = PsubE ./ RnormMax;
PobjEstar = PsubBsFun(PsubEstar,PsubBend,PsubMin);
ContkDA1scKrgEqHi = plot_contour(Felec, Pobj, Tsc5, PobjE(1), FscHMhi(1), 'off', 'ContkDA1scKrgEqHi', title_str1, ...
strcat({'Kda = '}, num2str(Kda(1))), graphs2files, FigDir);
ContkDA2scKrgEqHi = plot_contour(Felec, Pobj, Tsc6, PobjE(2), FscHMhi(2), 'off', 'ContkDA2scKrgEqHi', title_str2, ...
strcat({'Kda = '}, num2str(Kda(2))), graphs2files, FigDir);
bg_kDA1vskDA2scKrgEqHi = plot_bgStar(FscHM(1), FscHM(2), FscHMstar(1), FscHMstar(2),...
PobjE(1), PobjE(2), PobjEstar(1), PobjEstar(2),...
'off', 'kDA1vskDA2scKrgEqHi_bg', graphs2files, FigDir);
bg_root = 'bg_kDA1vskDA2scKrgEqHi';
quad_kDA1vskDA2scKrgEqHi = quad_subplot(ContkDA1scKrgEqHi, ContkDA2scKrgEqHi, bg_kDA1vskDA2scKrgEqHi, ...
'quad_kDA1vskDA2scKrgEqHi', bg_root, graphs2files, FigDir);
if show_graphics
quad_kDA1vskDA2scKrgEqHi.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "kDA1vskDA2_scKrgEqHi_quad", ...
"Effect of dopamine-transporter blockade on the reward mountain");
This section demontrates that the inconsistency between the simulated and observed results is seen over a broad range of values of the parameter that scales the MFB drive on the dopamine neurons.
The problem with the series-circuit model is a fundamental one. It predicts shifts of the mountain along the pulse-frequency axis in response to perturbation of dopaminergic neurotransmission. In contrast, systematic, consistent shifts along the pulse-frequency axis are not seen in eICSS studies under the influence of the dopamine transporter blocker, GBR-12909 (Hernandez et al. 2012); the dopamine, norepinephrine, and serotonin blocker, cocaine (Hernandez et al., 2010); the D2, D3 and 5HT7 receptor blocker, pimozide (Trujillo-Pisanty et al., 2014); or the cannabinoid CB-1 blocker, AM-251 (Trujillo-Pisanty et al., 2011). (AM-251 inhibits dopamine release and attenuates the stimulation-induced increase in dopamine tone (Trujillo-Pisanty et al., 2011)) . Thus, the series-circuit model fails to account for the eICSS data.
clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
The convergence model
tic;
close all;
The failure of the series-circuit model to account readily for the differential movement of the reward mountain in the eICSS and oICSS studies motivates the search for an alternative. We investigate a new model here. The challenge is to account for the stability of the reward mountain along the pulse-frequency axis under dopaminergic challenge in the eICSS data, the observed displacement along the pulse-frequency axis in the oICSS data, and the observed displacement along the price axis in both datasets.
if show_graphics
show_imported_graphic('eICSS+oICSS_convergence_noMasks_v11.png',30,ImpFigDir)
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "Convergence_model", ...
"A model implementing converging pathways subserving oICSS and eICSS");
In this model, the circuitry underlying eICSS of the medial forebrain bundle and oICSS of midbrain dopamine neurons includes parallel stages that are linked at two levels. The parallel limbs subserve oICSS and eICSS between the directly stimulated neurons and the scaled output of the reward-growth functions. The upstream link between the two limbs relays to midbrain dopamine neurons input from neurons activated by electrical stimulation of the medial forebrain bundle. The downsteam link combines the outputs of the two limbs so as to produce the signal representing the benefit of the experimenter-controlled reward.
The input from the MFB to the midbrain dopamine neurons
It has long been known that electrical stimulation of the medial forebrain bundle provides trans-synaptic input to midbrain dopamine neurons (Maeda & Mogenson, 1980) and drives phasic dopamine release from their ventral striatal terminals (Gratton, Hoffer & Gerhardt, 1988; Yavich & Tiihonen, 2000a, 2000b; Wightman & Robinson, 2002; Yavitch & Tanila, 2007). Recently, Cossette, Conover & Shizgal (2016) demonstrated differences between the frequency-following characteristics of the MFB input to midbrain dopamine neurons that project to the medial shell of the nucleus accumbens and the frequency-following characteristics of the neurons subserving eICSS of the MFB. Whereas the neurons subserving the rewarding effect maintain high-fidelity frequency following up to ~360 pulses (Solomon et al., 2015), the input to the nucleus-accumbens projecting dopamine neurons fails to follow pulse frequencies greater than 130 pulses . To accommodate this difference, a second low-pass filter is inserted below the link between the medial forebrain bundle fibers and the spike counter to the left of the dopamine neurons in
(Solomon et al., 2015), the input to the nucleus-accumbens projecting dopamine neurons fails to follow pulse frequencies greater than 130 pulses . To accommodate this difference, a second low-pass filter is inserted below the link between the medial forebrain bundle fibers and the spike counter to the left of the dopamine neurons indisp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=="Convergence_model")))}));
At pulse frequencies below 130 pulses , the amplitude of dopamine transients recorded in the nucleus accumbens in response to electrical stimulation of the medial forebrain bundle grows as function of both current and pulse frequency (Cossette, Conover & Shizgal, 2016). Increases in current can compensate for decreases in pulse frequency so as to hold constant the amplitude of the transient (and vice-versa). This suggests that with train duration held constant, the amplitude of the transient depends on the aggregate rate of firing in the medial forebrain bundle fibers that drive the activation of the dopamine neurons. That is the reason for the second spike counter.
, the amplitude of dopamine transients recorded in the nucleus accumbens in response to electrical stimulation of the medial forebrain bundle grows as function of both current and pulse frequency (Cossette, Conover & Shizgal, 2016). Increases in current can compensate for decreases in pulse frequency so as to hold constant the amplitude of the transient (and vice-versa). This suggests that with train duration held constant, the amplitude of the transient depends on the aggregate rate of firing in the medial forebrain bundle fibers that drive the activation of the dopamine neurons. That is the reason for the second spike counter. As in the case of the series-circuit model, the trans-synaptic drive on the dopamine neurons is represented in terms of the optical pulse frequency that provides equivalent dopaminergic activation.
Summation between the outputs of the two parallel limbs
The convergence model retains the idea that there is a final common path for neural signals that encode the predicted benefits of reward procurement. The two parallel circuit limbs, one subserving eICSS and the other oICSS, converge on this final common path. A simple way to implement this convergence is to add the outputs of the two limbs. This proposal faces a seemingly daunting challenge: electrical stimulation of the MFB activates midbrain dopamine neurons. If so, one would expect such stimulation to drive signaling in both of the hypothesized converging pathways. Wouldn't this produce at least some displacement of the reward mountain along the pulse-frequency axis in response to dopamine-transporter blockade? We show below that this is not necessarily the case. Indeed, the simulations show that given reasonable assumptions and values drawn from the current data, a convergence model can replicate the eICSS findings.
The reward-intensity signal at the output of the eICSS limb of the circuit is computed as above, by applying
disp(string({strcat({'Function '},num2str(fun_tab.Number(fun_tab.Name=='fRbsrFull')))}));
The reward-intensity signal at the output of the oICSS limb of the circuit is computed in the same manner as when the input consists of optical stimulation pulses. However, additional steps are required to compute the output of this limb during eICSS, when electrically excited medial forebrain bundle neurons provide the input to the dopamine neurons. The electrical pulse frequency is translated into the input to the dopamine neurons (expressed in terms of the equivalent optical pulse frequency) by the function, FmfbDAdriveFun.
% [fun_num, fun_tab] = add_fun(fun_num, fun_tab, ...
% "FmfbDAdriveFun", "FbendMFBda, Felec, FroMFBda, Kmfb, NnMFBs",...
% "Function that translates MFB drive on the DA neurons into an equivalent optical pulse frequency");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, ...
"FFmfbDAbend", "Felec, FFmfbDAmax, FoptEquivVec, FmfbDAro",...
"Function that translates MFB drive on the DA neurons into an equivalent optical pulse frequency");
We now simulate the output of the model in response to electrical stimulation of the medial forebrain bundle.
% Reward intensity produced by the upper limb
C = 0.473; % median from Sonnenschein et al., 2003
Delec = 0.5;
logFelecBend = 1.3222; % from Solomon et al., 2015
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587; % from Solomon et al., 2015
FelecRO = 10^logFelecRO;
NnMFB = 126;
RhoPi = 5000;
FPmax = 1000;
FhmUpper = FpulseHMfun(C, Delec, FelecBend, FPmax, FelecRO, NnMFB, RhoPi)
FhmUpper = 77.2222
numF = 121;
numParamVals = 6;
Felec = logspace(0,3,numF);
FelecMat = repmat(Felec,numParamVals,1);
gElec = 5;
KrgUpper = 1;
Rupper = fRbsr(Felec, FelecBend, FhmUpper, FelecRO, gElec, KrgUpper);
The position-parameter value for the upper limb (77 pulses ) is typical for the eICSS studies entailing measurement of the reward mountain.
) is typical for the eICSS studies entailing measurement of the reward mountain.disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='RGfunsElecFhm')))}));
TitleStrSemi = 'Upper-limb reward intensity';
pnam = "Fhm";
fnam = "Fhm";
% The data to be plotted must be in columns. Thus, FelecMat and RelecMat are transposed.
RG_MFBda_upper_semilog = plot_RG(Felec',Rupper',pnam,FhmUpper,fnam,TitleStrSemi,'lin', 'upper');
axh = findall(RG_MFBda_upper_semilog,'Type','Axes');
axh.XLabel.String = "Electrical pulse frequency";
axh.YLabel.String = "Upper-limb reward intensity";
if show_graphics
RG_MFBda_upper_semilog.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "RG_upper_semilog", ...
"Growth of upper-limb reward intensity");
Reward-intensity growth in the lower limb depends on the strength of the MFB drive,
% Reward intensity produced by the lower limb
% First compute the drive on the DA neurons in terms of the equivalent optical pulse frequency
FmfbDAbend = 40; % Make this gradual to match data in Cossette et al., 2016
FmfbDAro = 125; % Frequency following as described by Cossette et al., 2016
FFmfbDAmax = FilterFun(FPmax,FmfbDAbend, FmfbDAro);
FoptEquivVec = [5;10;20;40;80;160];
FoptEquivMat = repmat(FoptEquivVec, 1, numF);
FmfbDAdriveMat = FmfbDAdriveFun(FmfbDAbend, Felec, FFmfbDAmax, FoptEquivVec, FmfbDAro);
which is shown here:
% Columnar data are required to plot multiple lines. Thus the input matrices have been transposed.
FF_graphMFBda = plot_freqFoll(FelecMat',FmfbDAdriveMat','MFB_DA',1,1000,1,250);
% modify graph
lgnd = legend(num2str(FoptEquivVec),'Location','best');
lgnd.Title.String = "FoptEquiv";
axh = findall(FF_graphMFBda,'Type','Axes');
axh.XLabel.String = "Electrical pulse frequency";
axh.YLabel.String = "Equivalent optical pulse frequency";
if show_graphics
FF_graphMFBda.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "FreqFolMFBda", ...
"Electrically induced firing frequency in MFB input to the dopamine neurons");
D = Delec;
FdaBend = 20;
FdaRO = 50;
gDA = 5;
KrgLower = 1;
NnDA = 100; % Same value as in the series-circuit section
RhoPiDA = 1572; % Same value as in the series-circuit section for D = 0.5
FPmax = 1000;
FdaHM = FpulseHMfun(C, D, FdaBend, FPmax, FdaRO, NnDA, RhoPiDA)
FdaHM = 37.6342
The FdaHM parameter positions the reward-growth curve for the dopamine neurons along the axis representing the dopamine firing frequency:
logFopt = 0:0.025:2.4;
Fopt = 10.^logFopt;
Rlower = fRbsr(Fopt, FdaBend, FdaHM, FdaRO, gDA, KrgLower);
TitleStrSemi = 'Lower limb reward-intensity';
pnam = "FdaHM";
fnam = "FdaHM";
% The data to be plotted must be in columns. Thus, FelecMat and RelecMat are transposed.
RG_DA_LowerEquiv_semilog = plot_RG(Fopt',Rlower',pnam,FdaHM,fnam,TitleStrSemi,'lin', 'LowerEquiv');
RG_DA_LowerEquiv_semilog = modify_2D_graph(RG_DA_LowerEquiv_semilog, 'Axes', 'XLim', [1,250], ...
'RG_FdaHM_LowerEquiv_semilog', graphs2files, FigDir);
axh = findall(RG_DA_LowerEquiv_semilog,'Type','Axes');
axh.XLabel.String = "Equivalent optical pulse frequency";
axh.YLabel.String = "Lower-limb reward intensity";
if show_graphics
RG_DA_LowerEquiv_semilog.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "RG_DA_LowerEquiv_semilog", ...
"Lower-limb reward intensity as a function of equivalent optical pulse frequencies");
In the following graph, we see how reward-intensity grows in the lower limb as a function of the electrical pulse frequency.
Rlower = fRbsr(FmfbDAdriveMat, FdaBend, FdaHM, FdaRO, gDA, KrgLower);
TitleStrSemi = 'Lower-limb reward intensity';
pnam = "optDAeq";
fnam = "optDAeq";
% The data to be plotted must be in columns. Thus, FelecMat and RelecMat are transposed.
RG_MFBda_lower_semilog = plot_RG(FelecMat',Rlower',pnam,FoptEquivVec,fnam,TitleStrSemi,'lin', 'lower');
RG_MFBda_lower_semilog = modify_2D_graph(RG_MFBda_lower_semilog, 'Axes', 'XLim', [5,1000], ...
'RG_optDAeq_lower_semilog', graphs2files, FigDir);
axh = findall(RG_MFBda_lower_semilog,'Type','Axes');
axh.XLabel.String = "Electrical pulse frequency";
axh.YLabel.String = "Lower-limb reward intensity";
if show_graphics
RG_MFBda_lower_semilog.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "RG_lower_semilog", ...
"Growth of lower-limb reward intensity at multiple values of Kmfb");
The FmfbDAdrive function expresses the MFB drive on the dopamine neurons in units of equivalent optical pulses, i.e., an FmfbDAdrive of 25 means that the effect of the trans-synaptic input from the MFB is equivalent to that produced by 25 pulses of direct optical stimulation.
of direct optical stimulation.The graph above shows that strong input is required to generate a substantial response from the lower limb. The reason for this is that the value is fairly close to the assumed frequency-following limit of the dopamine cells. The parameter values used in this simulation yield a value close to the observed median in the current study when the train duration was 1 s. An adjustment is then made to predict model output for a train duration of 0.5 s, which was used in the eICSS work. This calculation are described above (see: The location parameter of the reward-growth function for BSR).
value is fairly close to the assumed frequency-following limit of the dopamine cells. The parameter values used in this simulation yield a value close to the observed median in the current study when the train duration was 1 s. An adjustment is then made to predict model output for a train duration of 0.5 s, which was used in the eICSS work. This calculation are described above (see: The location parameter of the reward-growth function for BSR).Growth curves for the summated reward intensity are shown here, at different values of the strength of the MFB drive on the dopamine neurons. The equivalent optical pulse frequencies (FoptEquiv) are the maximum values that the simulated MFB drive can achieve.
Rupper = repmat(Rupper,numParamVals,1); % add rows to Rupper so that it is the same size as Rlower
Rsum = Rupper + Rlower;
TitleStrSemi = 'Summated reward intensity';
pnam = "FoptEquiv";
fnam = "MFBsum";
% The data to be plotted must be in columns. Thus, FelecMat and RelecMat are transposed.
RG_MFBda_semilog = plot_RG(FelecMat',Rsum',pnam,FoptEquivVec,fnam,TitleStrSemi,'lin');
RG_MFBda_semilog = modify_2D_graph(RG_MFBda_semilog, 'Axes', 'XLim', [5,1000], ...
'RG_MFBda_semilog', graphs2files, FigDir);
axh = findall(RG_MFBda_semilog, 'Type', 'Axes');
axh.XLabel.String = "Electrical pulse frequency";
axh.YLabel.String = "Upper + lower reward intensity";
if show_graphics
RG_MFBda_semilog.Visible = 'on';
end
The reduced slope of the early portion of the leftmost curve is due to the fact that reward intensity in the lower limb begins rising at lower MFB pulse frequencies than reward intensity in the upper limb when the drive on the dopamine neurons is highest.
% Section break here to force display of RG_MFBda_semilog
close all;
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "RG_Kmfb_semilog", ...
"Growth of summated reward intensity at at multiple values of Kmfb");
To prepare for simulation of the mountain generated by the convergence model, we first derive the electrical pulse frequencies that produce half-maximal, summated reward intensity for each value of FoptEquiv
% Find FsumHM by means of interpolation
numF = 121;
Felec = logspace(0,3,numF);
[FsumHM, RsumMax] = find_FhmSum(Felec, Rsum);
We generate a mountain surface, scaling the maximum MFB drive on the dopamine neurons to be equivalent to an optical pulse frequency of 80 pulses , a maximal or near-maximal value in the current study.
, a maximal or near-maximal value in the current study. disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='RG_Kmfb_semilog')))}));
FoptEquiv = 80;
RupperFoptEquiv = Rupper(FoptEquivVec==FoptEquiv,:); % Selects the appropriate row from the matrix
RupperMax = max(RupperFoptEquiv);
RlowerFoptEquiv = Rlower(FoptEquivVec==FoptEquiv,:); % Selects the appropriate row from the matrix
RlowerMax = max(RlowerFoptEquiv);
Rsum = Rsum(FoptEquivVec==FoptEquiv,:);
FsumHM = FsumHM(FoptEquivVec==FoptEquiv);
RsumMax = RsumMax(FoptEquivVec==FoptEquiv);
FPmax = 1000;
RupperNormMax = fRbsrNorm(FPmax, FelecBend, FhmUpper, FelecRO, gElec);
RlowerNormMax = fRbsrNorm(FoptEquiv, FdaBend, FdaHM, FdaRO, gDA);
RupperRatio = RupperMax / (RlowerMax + RupperMax)
RupperRatio = 0.5671
RlowerRatio = RlowerMax / (RlowerMax + RupperMax)
RlowerRatio = 0.4329
RsumNormMax = (RupperRatio * RupperNormMax) + (RlowerRatio * RlowerNormMax); % Weighted average
dotPhiObj = 1;
Kaa = 1;
Kec = 1;
Krg = 1;
dotRaa = 0.1;
pObj = 1;
PsubBend = 0.5;
PsubMin = 1.82;
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RsumNormMax)
PobjE = 8.9716
a = 3;
numP = numF;
Pobj = logspace(0,3,numP); % row variable
Tsum = TAsumFun(a, Pobj, PobjE, PsubBend, PsubMin, Rsum', RsumNormMax);
% n.b., numel(Rsumxx) = numel(Pobj). Rsumxx has been transposed. Thus, Tsumxx is a square matrix.
MTNsum = plot_MTN(Felec, Pobj, Tsum, 'off', strcat('MTNsum',num2str(FoptEquiv)), 'reward mountain', ...
graphs2files, FigDir);
if show_graphics
MTN.Visible = 'on';
end
We now simulate the effect of dopamine-transporter blockade.
clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
tic;
close all;
The effect of dopamine transporter blockade in the convergence model
The reward-intensity signal produced by the upper and lower limbs of the circuit shown in
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=="Convergence_model")))}));
will first be computed separately and then combined additively, with two different values of , one representing the vehicle condition, and the second representing the condition in which the dopamine-transporter blocker has been administered .
, one representing the vehicle condition, and the second representing the condition in which the dopamine-transporter blocker has been administered .% Reward intensity produced by the upper limb
C = 0.473; % median from Sonnenschein et al., 2003
Delec = 0.5;
logFelecBend = 1.3222; % from Solomon et al., 2015
FelecBend = 10^logFelecBend;
logFelecRO = 2.5587; % from Solomon et al., 2015
FelecRO = 10^logFelecRO;
NnMFB = 126;
RhoPi = 5000;
FPmax = 1000;
FhmUpper = FpulseHMfun(C, Delec, FelecBend, FPmax, FelecRO, NnMFB, RhoPi)
FhmUpper = 77.2222
numF = 121;
numParamVals = 2; % for two values of Kda
Felec = logspace(0,3,numF);
FelecMat = repmat(Felec,numParamVals,1);
gElec = 5;
KrgUpper = 1;
Rupper = fRbsr(Felec, FelecBend, FhmUpper, FelecRO, gElec, KrgUpper);
% Reward intensity produced by the lower limb
% First compute the drive on the DA neurons in terms of the equivalent optical pulse frequency
FmfbDAbend = 40;
FmfbDAro = 125;
FFmfbDAmax = FilterFun(FPmax,FmfbDAbend, FmfbDAro);
FoptEquiv = 80; % maximal or near-maximal value in the current study
FmfbDAdrive = FmfbDAdriveFun(FmfbDAbend, Felec, FFmfbDAmax, FoptEquiv, FmfbDAro);
% Calculate FdaHM for the vehicle and drug conditions
C = 0.473; % median from Sonnenschein et al., 2003
D = Delec;
FdaBend = 20;
FdaRO = 50;
Kda = [1,10^0.15]; % Median for current oICSS study: 10^0.1446
gDA = 5; % Median in vehicle condition of current oICSS study: 4.61
KrgLower = 1;
NnDA = 100; % Make sure to use the same value as in the series-circuit section
RhoPiDA = 1572; % Same value as in the series-circuit section for D=0.5
FPmax = 1000;
FdaHMkDA = FpulseHMfun(C, D, FdaBend, FPmax, FdaRO, NnDA, RhoPiDA, Kda) % result is a two-element vector
% Compute reward growth in the lower limb as a function of the equivalent optical pulse frequency
numFOpt = 97; % 0.025 log10 spacing from 0 to 2.4
numParamVals = 2; % 2 values of Kda
Fopt = logspace(0,2.4,numFOpt);
FoptMat = repmat(Fopt,numParamVals,1); % number of rows = number of Kda vals
KdaMat = repmat(Kda',1,numFOpt); % make matrix the same size as FoptMat
FdaHMmat = repmat(FdaHMkDA',1,numFOpt); % make matrix the same size as FoptMat
Rlower = fRbsr(FoptMat, FdaBend, FdaHMmat, FdaRO, gDA, KrgLower);
TitleStrSemi = 'Lower limb vs Fopt equivalent';
pnam = "Kda";
fnam = "Kda";
% The data to be plotted must be in columns. Thus, FoptMat and Rlower are transposed.
RG_MFB_DA_FoptEquiv_Kda = plot_RG(FoptMat',Rlower',pnam,Kda,fnam,TitleStrSemi,'lin', 'FoptEquiv');
RG_MFB_DA_FoptEquiv_Kda = modify_2D_graph(RG_MFB_DA_FoptEquiv_Kda, 'Axes', 'XLim', [1,250], ...
'RG_Kda_FoptEquiv_semilog', graphs2files, FigDir);
axh = findall(RG_MFB_DA_FoptEquiv_Kda, 'Type','Axes');
axh.XLabel.String = "Equivalent optical pulse frequency";
axh.YLabel.String = "Lower-limb reward intensity";
if show_graphics
RG_MFB_DA_FoptEquiv_Kda.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "RG_MFB_DA_FoptEquiv_Kda", ...
"Reward growth in lower limb as a function of equivalent optical pulse frequency");
The previous graph shows lower-limb reward intensity as a function of the equivalent optical pulse frequency in the presence and absence of dopamine-transporter blockade. The next graph will plot lower-limb reward intensity as a function of the electrical MFB pulse frequency, again in the presence and absence of dopamine transporter blockade.
KdaMat = repmat(Kda',1,numF); % make matrix the same size as FelecMat
FdaHMmat = repmat(FdaHMkDA',1,numF); % make matrix the same size as FelecMat
Rlower = fRbsr(FmfbDAdrive, FdaBend, FdaHMmat, FdaRO, gDA, KrgLower);
TitleStrSemi = 'Lower limb vs Felec';
pnam = "Kda";
fnam = "Kda";
% The data to be plotted must be in columns. Thus, FelecMat and RelecMat are transposed.
RG_MFB_DA_Felec_Kda = plot_RG(FelecMat',Rlower',pnam,Kda,fnam,TitleStrSemi,'lin', 'Felec');
RG_MFB_DA_Felec_Kda = modify_2D_graph(RG_MFB_DA_Felec_Kda, 'Axes', 'XLim', [5,500], ...
'RG_Kda_Felec_semilog', graphs2files, FigDir);
axh = findall(RG_MFB_DA_Felec_Kda, 'Type','Axes');
axh.XLabel.String = "Electrical pulse frequency";
axh.YLabel.String = "Lower-limb reward intensity";
if show_graphics
RG_MFB_DA_Felec_Kda.Visible = 'on';
end
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, "RG_MFB_DA_Felec_Kda", ...
"Reward growth in lower limb as a function of electrical pulse frequency");
In the convergence model, the reward-intensity signal in the upper limb is not changed by dopamine-transporter blockade. Thus, we can simply add the reward-intensity signal already computed to the results for the lower-limb reward-intensity signal to obtain the summated output of the two limbs. This corresponds to the 'benefit' signal at the right of
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=="Convergence_model")))}));
Rupper = repmat(Rupper,numParamVals,1); % add rows to Rupper so that it is the same size as Rlower
Rsum = Rupper + Rlower;
The following graph shows the summated reward-intensity signal in response to MFB drive equivalent to an optical pulse frequency of
FoptEquiv
FoptEquiv = 80
pulses
[FsumHM, RsumMax] = find_FhmSum(FelecMat, Rsum)
% plot Rupper, Rlower, and Rxor for the vehicle and drug conditions
% Three curves will be plotted per graph; the matrices must be dimensioned accordingly
FelecMatPlot = repmat(Felec, 3, 1);
RmatVeh = [Rupper(1,:);Rlower(1,:);Rsum(1,:)];
RmatDrg = [Rupper(2,:);Rlower(2,:);Rsum(2,:)];
pnam = "1:up 2:lo 3:sum";
pVec = [1,2,3];
fnam = strcat("Sum",num2str(FoptEquiv),"Veh");
TitleStr = strcat("FoptEquiv=",num2str(FoptEquiv),"; vehicle");
RGsum_Rveh = plot_RG(FelecMatPlot',RmatVeh',pnam,pVec,fnam,TitleStr,'lin');
axh = findall(RGsum_Rveh, 'Type', 'Axes');
axh.XLabel.String = 'Electrical pulse frequency';
lh = findall(RGsum_Rveh, 'Type', 'Line');
lcolor = {'m','g','c'};
for j=1:3
lh(j).Color = lcolor{j};
end
if show_graphics
RGsum_Rveh.Visible = 'on';
end
fig_nam = strcat("RGsum_Veh_FoptEquiv",num2str(FoptEquiv));
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, fig_nam, "Reward summation in the convergence model");
fnam = strcat("sum",num2str(FoptEquiv),"Drg");
TitleStr = strcat("FoptEquiv=",num2str(FoptEquiv),"; drug");
RGsum_Rdrg = plot_RG(FelecMatPlot',RmatDrg',pnam,pVec,fnam,TitleStr,'lin');
axh = findall(RGsum_Rdrg, 'Type', 'Axes');
axh.XLabel.String = 'Electrical pulse frequency';
lh = findall(RGsum_Rdrg, 'Type', 'Line');
lcolor = {'m','g','c'};
for j=1:3
lh(j).Color = lcolor{j};
lh(j).LineStyle = '--';
end
if show_graphics
RG_RDrg.Visible = 'on';
end
fig_nam = strcat("RGsum_Drg_FoptEquiv",num2str(FoptEquiv));
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, fig_nam, "Reward summation in the convergence model");
shg
pnam = "Kda";
pVec = [Kda(1),Kda(2)];
fnam = strcat("Sum",num2str(FoptEquiv),"VehDrg");
TitleStr = strcat("FoptEquiv=",num2str(FoptEquiv),"; sum");
RGsum_Rvehdrg = plot_RG(FelecMat',Rsum',pnam,pVec,fnam,TitleStr,'lin');
axh = findall(RGsum_Rvehdrg, 'Type', 'Axes');
axh.XLabel.String = 'Electrical pulse frequency';
lh = findall(RGsum_Rvehdrg, 'Type', 'Line');
for j=1:2
lh(j).Color = 'm';
end
lh(1).LineStyle = '--';
lgndh = findall(RGsum_Rvehdrg,'Type','Legend');
lgndh.Title.String = 'condition';
lgndh.String = {'vehicle','drug'};
if show_graphics
RGsum_Rvehdrg.Visible = 'on';
end
fig_nam = strcat("RGsum_VehDrg_FoptEquiv",num2str(FoptEquiv));
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, fig_nam, "Drug modulation of reward summation in the convergence model");
At this high value of MFB drive, a small leftward curve shift is produced.
With the reward-intensity signal in hand, we can proceed to generate the mountains predicted in the presence and absence of dopamine-transporter blockade.
close all;
FPmax = 1000;
FhmUpperVec = repmat(FhmUpper,1,2); % Upper limb output is independent of Kda
RupperNormMax = fRbsrNorm(FPmax, FelecBend, FhmUpperVec, FelecRO, gElec)';
RlowerNormMax = fRbsrNorm(FoptEquiv, FdaBend, FdaHMkDA, FdaRO, gDA)';
RupperMax = max(Rupper, [], 2);
RlowerMax = max(Rlower, [], 2);
RupperRatio = RupperMax ./ (RlowerMax + RupperMax)
RlowerRatio = RlowerMax ./ (RlowerMax + RupperMax)
RsumNormMax = (RupperRatio .* RupperNormMax) + (RlowerRatio .* RlowerNormMax); % Weighted average
dotPhiObj = 1;
Kaa = 1;
Kec = 1;
Krg = 1;
dotRaa = 0.1;
pObj = 1;
PsubBend = 0.5;
PsubMin = 1.82;
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RsumNormMax)
a = 3;
numP = numF;
Pobj = logspace(0,3,numP); % row variable
TsumKda1 = TAsumFun(a, Pobj, PobjE(1), PsubBend, PsubMin, Rsum(1,:)', RsumNormMax(1));
TsumKda2 = TAsumFun(a, Pobj, PobjE(2), PsubBend, PsubMin, Rsum(2,:)', RsumNormMax(2));
% n.b., numel(Rsum(x)) = numel(Pobj). Rsum(x) has been transposed. Thus, TsumKda(x) is a square matrix.
% Plot the individual mountains for the two values of Kda
title_str1 = strcat({'Kda = '}, sprintf('%2.2f', Kda(1)));
MTNsumKda1 = plot_MTN(Felec, Pobj, TsumKda1, 'off', 'MTNsumKda1', title_str1, ...
graphs2files, FigDir);
title_str2 = strcat({'Kda = '}, sprintf('%2.2f', Kda(2)));
MTNsumKda2 = plot_MTN(Felec, Pobj, TsumKda2, 'off', 'MTNsumKda1', title_str2, ...
graphs2files, FigDir);
% Plot the pair of mountains for the two values of Kda
dual_sum_plot = dual_subplot(MTNsumKda1, MTNsumKda2, 'MTNsum_Kda1_Kda2',...
graphs2files,FigDir);
if show_graphics
dual_sum_plot.Visible = 'on';
end
% Estimate the weighted average of the firing frequencies that produce half-maximal reward intensity in the
% upper and lower limbs
FhmUpperStar = FilterFun(FhmUpper,FelecBend,FelecRO);
FhmLowerDrive = FmfbDAdriveFun(FmfbDAbend, FsumHM, FFmfbDAmax, FoptEquiv, FmfbDAro);
FFhmLower = FilterFun(FhmLowerDrive,FdaBend, FdaRO);
FsumHMstar = FsumHM .* ...
((RupperRatio .* (FhmUpperStar ./ FsumHM)) + (RlowerRatio .* (FFhmLower ./ FhmLowerDrive)));
% Weighted average of the frequency-following fidelity in the upper and lower limbs
PsubE = PsubFun(PobjE,PsubBend,PsubMin);
PsubEstar = PsubE ./ RsumNormMax;
PobjEstar = PsubBsFun(PsubEstar,PsubBend,PsubMin);
ContkDA1sum = plot_contour(Felec, Pobj, TsumKda1, PobjE(1), FsumHM(1), 'off', 'ContkDA1sum', title_str1, ...
strcat({'Kda = '}, num2str(Kda(1))), graphs2files, FigDir);
ContkDA2sum = plot_contour(Felec, Pobj, TsumKda2, PobjE(2), FsumHM(2), 'off', 'ContkDA2sum', title_str2, ...
strcat({'Kda = '}, num2str(Kda(2))), graphs2files, FigDir);
bg_kDA1vskDA2sum = plot_bgStar(FsumHM(1), FsumHM(2), FsumHMstar(1), FsumHMstar(2),...
PobjE(1), PobjE(2), PobjEstar(1), PobjEstar(2),...
'off', 'kDA1vskDA2sum_bg', graphs2files, FigDir, -0.3, 0.3);
% Set a scale for the y-axis of this bar graph that ressembles the scales for the others
bg_root = 'bg_kDA1vskDA2sum';
quad_kDA1vskDA2sum = quad_subplot(ContkDA1sum, ContkDA2sum, bg_kDA1vskDA2sum, 'quad_kDA1vskDA2sum', bg_root, ...
graphs2files, FigDir);
if show_graphics
quad_kDA1vskDA2sum.Visible = 'on';
end
fig_nam = strcat("quad_kDA1vskDA2sum",num2str(FoptEquiv));
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, ...
fig_nam, "Modest downward shift of the reward mountain induced by dopamine-transporter blockade");
% Insert section break to force display of the last graph
close all;
As expected from the small leftward displacement of the reward-intensity growth function by dopamine transporter blockade shown in
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=="RGsum_VehDrg_FoptEquiv80")))}));
the mountain is shifted modestly along the pulse-frequency axis. None of the eight subjects in the Hernandez et al. (2012) eICSS study showed such a shift.
Adding dopaminergic modulation of effort cost (i.e., a drug-induced reduction in these costs) boosts the shift along the price axis. The same effect could be produced by boosting the output of the reward-growth function(s) 
or decreasing the value of alternate activities.
or decreasing the value of alternate activities.With dopaminergic modulation of effort cost:
KeffMod = [1;10^0.15]; % Median shift in Pe in the current oICSS study: 10^0.1525
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RsumNormMax, KeffMod)
TsumEffModKda1 = TAsumFun(a, Pobj, PobjE(1), PsubBend, PsubMin, Rsum(1,:)', RsumNormMax(1));
TsumEffModKda2 = TAsumFun(a, Pobj, PobjE(2), PsubBend, PsubMin, Rsum(2,:)', RsumNormMax(2));
% n.b., numel(Rsum(x)) = numel(Pobj). Rsum(x) has been transposed. Thus, TsumKda(x) is a square matrix.
% Plot the individual mountains for the two values of Kda
title_str1 = strcat({'Kda = '}, sprintf('%2.2f', Kda(1)));
MTNsumEffModKda1 = plot_MTN(Felec, Pobj, TsumEffModKda1, 'off', 'MTNsumEffModKda1', title_str1, ...
graphs2files, FigDir);
title_str2 = strcat({'Kda = '}, sprintf('%2.2f', Kda(2)));
MTNsumEffModKda2 = plot_MTN(Felec, Pobj, TsumEffModKda2, 'off', 'MTNsumEffModKda2', title_str2, ...
graphs2files, FigDir);
% Plot the pair of mountains for the two values of Kda
dual_sum_plot = dual_subplot(MTNsumEffModKda1, MTNsumEffModKda2, 'MTNsumEffMod_Kda1_Kda2',...
graphs2files,FigDir);
if show_graphics
dual_sum_plot.Visible = 'on';
end
% Use value of FsumHMstar calculated above;
PsubE = PsubFun(PobjE,PsubBend,PsubMin);
PsubEstar = PsubE ./ RsumNormMax;
PobjEstar = PsubBsFun(PsubEstar,PsubBend,PsubMin);
ContkDA1sumEffMod = plot_contour(Felec, Pobj, TsumEffModKda1, PobjE(1), FsumHM(1), 'off', ...
'ContkDA1sumEffMod', title_str1, ...
strcat({'Kda = '}, num2str(Kda(1))), graphs2files, FigDir);
ContkDA2sumEffMod = plot_contour(Felec, Pobj, TsumEffModKda2, PobjE(2), FsumHM(2), 'off', ...
'ContkDA2sumEffMod', title_str2, ...
strcat({'Kda = '}, num2str(Kda(2))), graphs2files, FigDir);
bg_kDA1vskDA2sumEffMod = plot_bgStar(FsumHM(1), FsumHM(2), FsumHMstar(1), FsumHMstar(2),...
PobjE(1), PobjE(2), PobjEstar(1), PobjEstar(2),...
'off', 'kDA1vskDA2sumEffMod_bg', graphs2files, FigDir, -0.3, 0.3);
% Set a scale for the y-axis of this bar graph that ressembles the scales for the others
bg_root = 'bg_kDA1vskDA2sumEffMod';
quad_kDA1vskDA2sumEffMod = quad_subplot(ContkDA1sumEffMod, ContkDA2sumEffMod, bg_kDA1vskDA2sumEffMod, ...
'quad_kDA1vskDA2sumEffMod', bg_root, graphs2files, FigDir);
if show_graphics
quad_kDA1vskDA2sumEffMod.Visible = 'on';
end
fig_nam = strcat("quad_kDA1vskDA2sumEffMod",num2str(FoptEquiv));
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, ...
fig_nam, "Dual shifts of the reward mountain induced by dopamine-transporter blockade");
close all;
vars2save = who;
ws_name = fullfile(pwd, 'WS_convergence_vars80.mat');
save_ws;
The small shift along the pulse-frequency axis can be eliminated by reducing the MFB drive. We will illustrate this by means of reward-growth graphs.
% Rupper is unchanged - no need to regenerate it
% Regenerate Rlower using weaker MFB drive
FoptEquiv = 40; % Value was 80 in above simulations
FmfbDAdrive = FmfbDAdriveFun(FmfbDAbend, Felec, FFmfbDAmax, FoptEquiv, FmfbDAro);
FdaHMkDA = FpulseHMfun(C, D, FdaBend, FPmax, FdaRO, NnDA, RhoPiDA, Kda) % result is a two-element vector
KdaMat = repmat(Kda',1,numF);
FdaHMmat = repmat(FdaHMkDA',1,numF);
Rlower = fRbsr(FmfbDAdrive, FdaBend, FdaHMmat, FdaRO, gDA, KrgLower);
Rsum = Rupper + Rlower;
[FsumHM, RsumMax] = find_FhmSum(FelecMat, Rsum)
RmatVeh = [Rupper(1,:);Rlower(1,:);Rsum(1,:)];
RmatDrg = [Rupper(2,:);Rlower(2,:);Rsum(2,:)];
pnam = "1:up 2:lo 3:sum";
pVec = [1,2,3];
fnam = strcat("Sum",num2str(FoptEquiv),"Veh");
TitleStr = strcat("FoptEquiv=",num2str(FoptEquiv),"; vehicle");
RGsum_Rveh = plot_RG(FelecMatPlot',RmatVeh',pnam,pVec,fnam,TitleStr,'lin');
axh = findall(RGsum_Rveh, 'Type', 'Axes');
axh.XLabel.String = 'Electrical pulse frequency';
lh = findall(RGsum_Rveh, 'Type', 'Line');
lcolor = {'m','g','c'};
for j=1:3
lh(j).Color = lcolor{j};
end
if show_graphics
RGsum_Rveh.Visible = 'on';
end
fig_nam = strcat("RGsum_Veh_FoptEquiv",num2str(FoptEquiv));
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, fig_nam, "Reward summation in the convergence model");
pnam = "1:up 2:lo 3:sum";
pVec = [1,2,3];
fnam = strcat("Sum",num2str(FoptEquiv),"Drg");
TitleStr = strcat("FoptEquiv=",num2str(FoptEquiv),"; drug");
RGsum_Rdrg = plot_RG(FelecMatPlot',RmatDrg',pnam,pVec,fnam,TitleStr,'lin');
axh = findall(RGsum_Rdrg, 'Type', 'Axes');
axh.XLabel.String = 'Electrical pulse frequency';
lh = findall(RGsum_Rdrg, 'Type', 'Line');
lcolor = {'m','g','c'};
for j=1:3
lh(j).Color = lcolor{j};
lh(j).LineStyle = '--';
end
if show_graphics
RGsum_Rdrg.Visible = 'on';
end
fig_nam = strcat("RGsum_Drg_FoptEquiv",num2str(FoptEquiv));
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, fig_nam, "Reward summation in the convergence model");
pnam = "Kda";
pVec = [Kda(1),Kda(2)];
fnam = strcat("Sum",num2str(FoptEquiv),"VehDrg");
TitleStr = strcat("FoptEquiv=",num2str(FoptEquiv),"; sum");
RGsum_Rvehdrg = plot_RG(FelecMat',Rsum',pnam,pVec,fnam,TitleStr,'lin');
axh = findall(RGsum_Rvehdrg, 'Type', 'Axes');
axh.XLabel.String = 'Electrical pulse frequency';
lh = findall(RGsum_Rvehdrg, 'Type', 'Line');
lh(1).LineStyle = '--';
for j=1:2
lh(j).Color = 'm';
end
lgndh = findall(RGsum_Rvehdrg, 'Type', 'Legend');
lgndh.Title.String = "condition";
lgndh.String = {'vehicle','drug'};
if show_graphics
RGsum_Rvehdrg.Visible = 'on';
end
fig_nam = strcat("RGsum_VehDrg_FoptEquiv",num2str(FoptEquiv));
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, fig_nam, "Drug modulation of reward summation in the convergence model");
At this weaker value of MFB drive, the reward-growth curve is shifted to the right. At this value of (~median for the current study), the drug cannot shift the reward-growth curve far enough to the left to displace the summated curve leftwards. Thus, we learn from this that shifts along the pulse-frequency axis in the convergence model depend on the strength of the MFB drive.
(~median for the current study), the drug cannot shift the reward-growth curve far enough to the left to displace the summated curve leftwards. Thus, we learn from this that shifts along the pulse-frequency axis in the convergence model depend on the strength of the MFB drive.close all;
FPmax = 1000;
FhmUpperVec = repmat(FhmUpper,1,2); % Upper limb output is independent of Kda
RupperNormMax = fRbsrNorm(FPmax, FelecBend, FhmUpperVec, FelecRO, gElec)';
RlowerNormMax = fRbsrNorm(FoptEquiv, FdaBend, FdaHMkDA, FdaRO, gDA)';
RupperMax = max(Rupper, [], 2);
RlowerMax = max(Rlower, [], 2);
RupperRatio = RupperMax ./ (RlowerMax + RupperMax)
RlowerRatio = RlowerMax ./ (RlowerMax + RupperMax)
RsumNormMax = (RupperRatio .* RupperNormMax) + (RlowerRatio .* RlowerNormMax); % Weighted average
dotPhiObj = 1;
Kaa = 1;
Kec = 1;
Krg = 1;
dotRaa = 0.1;
pObj = 1;
PsubBend = 0.5;
PsubMin = 1.82;
PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, dotRaa, pObj, PsubBend, PsubMin, RsumNormMax, KeffMod)
a = 3;
numP = numF;
Pobj = logspace(0,3,numP); % row variable
TsumEffModLoFoeqKda1 = TAsumFun(a, Pobj, PobjE(1), PsubBend, PsubMin, Rsum(1,:)', RsumNormMax(1));
TsumEffModLoFoeqKda2 = TAsumFun(a, Pobj, PobjE(2), PsubBend, PsubMin, Rsum(2,:)', RsumNormMax(2));
% n.b., numel(Rsum(x)) = numel(Pobj). Rsum(x) has been transposed. Thus, TsumKda(x) is a square matrix.
% Plot the individual mountains for the two values of Kda
title_str1 = strcat({'Kda = '}, sprintf('%2.2f', Kda(1)));
MTNsumEffModLoFoeqKda1 = plot_MTN(Felec, Pobj, TsumEffModLoFoeqKda1, 'off', 'MTNsumEffModLoFoeqKda1', title_str1, ...
graphs2files, FigDir);
title_str2 = strcat({'Kda = '}, sprintf('%2.2f', Kda(2)));
MTNsumEffModKLoFoeqda2 = plot_MTN(Felec, Pobj, TsumEffModLoFoeqKda2, 'off', 'MTNsumEffModLoFoeqKda2', title_str2, ...
graphs2files, FigDir);
% Plot the pair of mountains for the two values of Kda
dual_sum_plot = dual_subplot(MTNsumEffModLoFoeqKda1, MTNsumEffModKLoFoeqda2, 'MTNsumEffModLoFoeq_Kda1_Kda2',...
graphs2files,FigDir);
if show_graphics
dual_sum_plot.Visible = 'on';
end
% Estimate the weighted average of the firing frequencies that produce half-maximal reward intensity in the
% upper and lower limbs
FhmUpperStar = FilterFun(FhmUpper,FelecBend,FelecRO);
FhmLowerDrive = FmfbDAdriveFun(FmfbDAbend, FsumHM, FFmfbDAmax, FoptEquiv, FmfbDAro);
FFhmLowerStar = FilterFun(FhmLowerDrive,FdaBend, FdaRO);
FsumHMstar = FsumHM .* ...
((RupperRatio .* (FhmUpperStar ./ FsumHM)) + (RlowerRatio .* (FFhmLowerStar ./ FhmLowerDrive)));
PsubE = PsubFun(PobjE,PsubBend,PsubMin);
PsubEstar = PsubE ./ RsumNormMax;
PobjEstar = PsubBsFun(PsubEstar,PsubBend,PsubMin);
ContkDA1sumEffMod = plot_contour(Felec, Pobj, TsumEffModLoFoeqKda1, PobjE(1), FsumHM(1), 'off', ...
'ContkDA1sumEffModLoFoeq', title_str1, ...
strcat({'Kda = '}, num2str(Kda(1))), graphs2files, FigDir);
ContkDA2sumEffMod = plot_contour(Felec, Pobj, TsumEffModLoFoeqKda2, PobjE(2), FsumHM(2), 'off', ...
'ContkDA2sumEffModLoFoeq', title_str2, ...
strcat({'Kda = '}, num2str(Kda(2))), graphs2files, FigDir);
bg_kDA1vskDA2sumEffMod = plot_bgStar(FsumHM(1), FsumHM(2), FsumHMstar(1), FsumHMstar(2),...
PobjE(1), PobjE(2), PobjEstar(1), PobjEstar(2),...
'off', 'kDA1vskDA2sumEffMod_bg', graphs2files, FigDir, -0.3, 0.3);
% Set a scale for the y-axis of this bar graph that ressembles the scales for the others
bg_root = 'bg_kDA1vskDA2sumEffMod';
quad_kDA1vskDA2sumEffMod = quad_subplot(ContkDA1sumEffMod, ContkDA2sumEffMod, bg_kDA1vskDA2sumEffMod, ...
'quad_kDA1vskDA2sumEffMod', bg_root, graphs2files, FigDir);
if show_graphics
quad_kDA1vskDA2sumEffMod.Visible = 'on';
end
fig_nam = strcat("quad_kDA1vskDA2sumEffMod",num2str(FoptEquiv));
[fig_num, fig_tab] = add_fig(fig_num, fig_tab, ...
fig_nam, "Dual shifts of the reward mountain induced by dopamine-transporter blockade");
close all;
vars2save = who;
ws_name = fullfile(pwd, 'WS_convergence_vars.mat');
save_ws;
clearvars('-except',keepVars{:})
keepVars = who; % Restore cell array containing names of variables to be retained
toc
tic;
Remarks concerning the convergence model
The results in
disp(string({strcat({'Figure '},num2str(fig_tab.Number(fig_tab.Name=='quad_kDA1vskDA2sumEffMod')))}));
above resemble those obtained in the eICSS study in which the reward mountain was measured under the influence of GBR-12901 (Hernandez et al., 2012). Thus, given moderate MFB drive and values compatible with those obtained in the present study (adjusted for a train duration of 0.5 s), the convergence model can generate outputs that match our earlier eICSS results.
values compatible with those obtained in the present study (adjusted for a train duration of 0.5 s), the convergence model can generate outputs that match our earlier eICSS results.The dependence of the output of the convergence model on experimental parameters is of interest. Although failure to observe shifts along the pulse-frequency axis was the most common result of the experiments in which the eICSS reward-mountain was measured under the influence of drugs that alter dopaminergic neurotransmission (27 of 32 cases reported by Hernandez et al., 2010, Trujillo-Pisanty et al., 2011, Hernandez et al., 2012, and Trujillo-Pisanty et al., 2014), it is not the only result. Reward mountains obtained from three subjects in the cocaine study (Hernandez et al., 2010) showed fairly substantial, reliable, shifts along the pulse-frequency axis. (The shift was marginally reliable in a fourth subject when tested initially but disappeared upon re-test.) Although no subject in the GBR-12909 (Hernandez et al., 2012) or pimozide (Trujillo-Pisanty et al., 2014) studies showed such shifts, one subject in the AM-251 study (Trujillo-Pisanty et al., 2011) did. Could variation in the strength of MFB drive on the dopamine neurons and in the individual values for the dopamine pathway explain these findings?
values for the dopamine pathway explain these findings?Please see the main body of the manuscript for more detailed discussion.
Sort tables and write to files
Formatted equations
if tabdir
writetable(sortrows(eqn_tab,2),...
fullfile(tabdir,strcat('eqn_tab_v',num2str(version),'.csv')),...
'WriteVariableNames', true);
end
Figures
if tabdir
writetable(sortrows(fig_tab,2),...
fullfile(tabdir,strcat('fig_tab_v',num2str(version),'.csv')),...
'WriteVariableNames', true);
end
Functions
if tabdir
writetable(sortrows(fun_tab,2),...
fullfile(tabdir,strcat('fun_tab_v',num2str(version),'.csv')),...
'WriteVariableNames', true);
end
Symbols
if tabdir
writetable(sortrows(sym_tab,2),...
fullfile(tabdir,strcat('sym_tab_v',num2str(version),'.csv')),...
'WriteVariableNames', true);
end
Save workspace
if saveWS
ws_name = fullfile(pwd, strcat('WS_GBR_eICSS_oICSS_v',num2str(version),'.mat'));
save_ws;
end
Last revised on:
datetime('now')
toc % elapsed time for this section
toc(t_start) % elapsed time for the entire script
Tools for this live script
function init_all
global eqn_num eqn_tab fig_num fig_tab fun_num fun_tab sym_num sym_tab
eqn_num = 0; % initialize counter
eqn_nam = ""; % initialize equation-name string array
eqn_desc = ""; % initialize equation-description string array;
eqn_tab = table(eqn_num, eqn_nam, eqn_desc,'VariableNames', {'Number'; 'Name'; 'Description'});
eqn_tab(1,:)=[]; % clear the table
disp("The equation table has been initialized.");
fig_num = 0; % initialize counter
fig_nam = ""; % initialize figure-name string array
fig_desc = ""; % initialize figure-description string array;
fig_tab = table(fig_num, fig_nam, fig_desc,'VariableNames', {'Number'; 'Name'; 'Description'});
fig_tab(1,:)=[]; % clear the table
disp("The figure table has been initialized.");
fun_num = 0; % initialize counter
fun_nam = ""; % initialize function-name string array
fun_args = ""; %initialize function-argument string array
fun_desc = ""; % initialize function-description string array;
fun_tab = table(fun_num, fun_nam, fun_args, fun_desc, ...
'VariableNames', {'Number'; 'Name'; 'Arguments'; 'Description'});
fun_tab(1,:)=[]; % clear the table
disp("The function table has been initialized.");
sym_num = 0; % initialize counter
sym_nam = ""; % initialize symbol-name string array
sym_desc = ""; % initialize symbol-description string array;
sym_tab = table(sym_num, sym_nam, sym_desc,'VariableNames', {'Number'; 'Name'; 'Description'});
sym_tab(1,:)=[]; % clear the table
disp("The symbol table has been initialized.");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "init_all", "", "initialize tables");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "add_eqn", "eqn_num, eqn_tab, EqnNam, EqnDesc",...
"add equation to table");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "add_fig", "fig_num, fig_tab, FigNam, FigDesc",...
"add figure to table");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "add_fun", "fun_num, fun_tab, FunNam, FunArgs, FunDesc",...
"add local function to table");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "add_sym", "sym_num, sym_tab, SymNam, SymDesc",...
"add symbol to table");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "FilterFun", "F, Fbend, Fro",...
"Frequency-following function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "FilterFunBS", "Fbend, Fro, Frate, FPmax",...
"Back-solution of frequency-following function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "LogistNormFun", "exponent, input, location",...
"Rising logistic function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "LogistNormBsFun", "exponent, location, output",...
"Back-solution of logistic function to return the input");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "LogistNormBsLocFun", "exponent, input, output",...
"Back-solution of logistic function to return the location parameter");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "PsubFun", "Pobj, PsubBend, PsubMin",...
"Subjective-price function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "PsubBsFun", "Psub, PsubBend, PsubMin",...
"Back-solution of subjective-price function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "ScalarDivFun", "dividend, divisor",...
"Scalar-division function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "ScalarDivBsFun", "divisor, quotient",...
"Back-solution of scalar-division function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "ScalarMultFun", "multiplicand, multiplier",...
"Scalar-multiplication function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "ScalarMultBsFun", "multiplier, product",...
"Back-solution of scalar-multiplication function");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "PsubEfun", "Ceff, Rmax, Ue",...
"Function to compute subjective equivalent of the objective-price location parameter");
[fun_num, fun_tab] = add_fun(fun_num, fun_tab, "PobjEfun", "Ceff, PsubBend, PsubMin, Rmax, Ue",...
"Function to compute objective-price location parameter");
end % end of init_all function
function [eqn_num, eqn_tab] = add_eqn(eqn_num, eqn_tab, EqnNam, EqnDesc)
if any(strcmp(EqnNam,eqn_tab.Name))
disp(strcat({'Equation '}, EqnNam, {' has already been entered.'}));
return;
else
eqn_num = eqn_num+1;
eqn_tab(eqn_num,:) = {eqn_num, ...
string(EqnNam), ...
string(EqnDesc)};
disp(strcat({'This is equation #'} , num2str(eqn_tab.Number(height(eqn_tab))), ...
{': ('}, eqn_tab.Name(height(eqn_tab)), {')'}));
end
end
function [fig_num, fig_tab] = add_fig(fig_num, fig_tab, FigNam, FigDesc)
if any(strcmp(FigNam,fig_tab.Name))
disp(strcat({'Figure '}, FigNam, {' has already been entered.'}));
return;
else
fig_num = fig_num+1;
fig_tab(fig_num,:) = {fig_num, ...
string(FigNam), ...
string(FigDesc)};
disp(strcat({'This is figure #'} , num2str(fig_tab.Number(height(fig_tab))), ...
{': ('}, fig_tab.Name(height(fig_tab)), {')'}));
end
end
function [fun_num, fun_tab] = add_fun(fun_num, fun_tab, FunNam, FunArgs, FunDesc)
if any(strcmp(FunNam,fun_tab.Name))
disp(strcat({'Function '}, FunNam, {' has already been entered.'}));
return;
else
fun_num = fun_num+1;
fun_tab(fun_num,:) = {fun_num, ...
string(FunNam), ...
string(FunArgs), ...
string(FunDesc)};
disp(strcat({'This is function #'} , num2str(fun_tab.Number(height(fun_tab))), ...
{': ('}, fun_tab.Name(height(fun_tab)), {')'}));
end
end
function [sym_num, sym_tab] = add_sym(sym_num, sym_tab, SymNam, SymDesc)
if any(strcmp(SymNam,sym_tab.Name))
disp(strcat({'Symbol '}, SymNam, {' has already been entered.'}));
return;
else
sym_num = sym_num+1;
sym_tab(sym_num,:) = {sym_num, ...
string(SymNam), ...
string(SymDesc)};
% disp(strcat({'This is symbol #'} , num2str(sym_tab.Number(height(sym_tab))), ...
% {': ('}, sym_tab.Name(height(sym_tab)), {')'}));
end
end
function FigDir = set_figdir(figdirnam, graphs2files)
global Figdir
FigDir = fullfile(pwd,figdirnam);
if graphs2files
if ~exist(FigDir, 'dir')
mkdir(FigDir);
end
end
end
function ImpFigDir = set_impfigdir(impfigdirnam)
ImpFigDir = fullfile(pwd,impfigdirnam);
if ~exist(ImpFigDir, 'dir')
disp(strcat({'Exiting because Imported-files directory '}, impfigdirnam, {' does not exist.'}));
return
end
end
function tabdir = set_tabdir(tabdirnam, tabs2files)
tabdir = fullfile(pwd,tabdirnam);
if tabs2files
if ~exist(tabdir, 'dir')
mkdir(tabdir);
end
end
end
Functional building blocks for the simulations
function Frate = FilterFun(F, Fbend, Fro) % Frequency-following function
Frate = Fbend.*(log(1+exp(Fro./Fbend))-log(1+exp((Fro-F)./Fbend)));
end
% see:
% Breton, Y.-A., Mullett, A., Conover, K., & Shizgal, P. (2013).
% Validation and extension of the reward-mountain model.
% Frontiers in Behavioral Neuroscience, 7, 125. https://doi.org/10.3389/fnbeh.2013.00125
% Solomon, R. B., Trujillo-Pisanty, I., Conover, K., & Shizgal, P. (2015).
% Psychophysical inference of frequency-following fidelity in the neural substrate for brain stimulation reward.
% Brain Research, 292, 327–341. https://doi.org/10.1016/j.bbr.2015.06.008
% Last revised by Peter Shizgal on 2019-09-05 09:57
function Fpulse = FilterFunBS(Fbend, Fro, FF, FPmax) % Back-solution of frequency-following function
% Estimate the maximum achievable firing rate, given FPmax, Fbend, and Fro
FFmax = FilterFun(FPmax, Fbend, Fro);
FpulseRI = Fro - Fbend.*log(exp(-FF/Fbend).*(exp(Fro/Fbend) + 1) - 1);
% FpulseRI is a complex number. When FF < FFmax, the imaginary part = zero. When Fpulse RI = FFmax, then Fpulse = Inf.
% We retain the real component when FF < FFmax & >= 0. returning NaN when FF > FFmax | FF < 0 and Inf when FF = FFmax.
sz = size(FF);
Fpulse = zeros(sz); % pre-allocate
neg = FF < 0;
below = (FF < FFmax & FF >= 0);
eq = FF == FFmax;
above = FF > FFmax;
Fpulse(neg) = NaN;
Fpulse(below) = real(FpulseRI(below)); % Complex portion, if it exists, equals zero
Fpulse(eq) = Inf;
Fpulse(above) = NaN;
end
% This back-solution can be verified from the following code snippet, which requires the
% Matlab Symbolic-Math Toolbox. It is commented out here to allow this Live Script to run
% on systems that lack the Symbolic-Math toolbox.
% syms FdaRate FdaBend FdaRO Fopt positive;
% simplify(solve(FdaRate == FdaBend * (log(1+exp(FdaRO/FdaBend)) - (log(1+exp((FdaRO-Fopt)/FdaBend)))),Fopt))
% Fopt = FdaRO - FdaBend*log(exp(-FdaRate/FdaBend)*(exp(FdaRO/FdaBend) + 1) - 1)
%%
% This function returns imaginary numbers. The imaginary portion equals zero provided that Frate (the firing rate)
% is less than its maximum achievable value. That maximum is estimated from the forward solution, using FPmax
% as the input.
%%
% see:
% Breton, Y.-A., Mullett, A., Conover, K., & Shizgal, P. (2013).
% Validation and extension of the reward-mountain model.
% Frontiers in Behavioral Neuroscience, 7, 125. https://doi.org/10.3389/fnbeh.2013.00125
% Solomon, R. B., Trujillo-Pisanty, I., Conover, K., & Shizgal, P. (2015).
% Psychophysical inference of frequency-following fidelity in the neural substrate for brain stimulation reward.
% Behavioural Brain Research, 292, 327–341. https://doi.org/10.1016/j.bbr.2015.06.008
function output = LogistNormFun(exponent, input, location)
output = ...
((input .^ exponent) ./...
((input .^ exponent) + (location .^ exponent)));
end
% The symbol for the logistic function in the flow diagrams is a rising S-shaped curve in a box.
%
% Revised by Peter Shizgal on 2019-10-29 21:25
% The argument called "output" must be scalar
function input = LogistNormBsFun(exponent, location, output)
% Logistic function back-solved to return the input
if (output >= 1) || (output <= 0)
disp("Illegal value. Call with 0 < output < 1");
return
end
input = location .* (output./(1-output)).^(1./exponent);
end
% Revised by Peter Shizgal on 2019-10-29 21:25
function location = LogistNormBsLocFun(exponent, input, output)
% The argument called "output" must be scalar
% Logistic function back-solved to return the location parameter
if output >= 1 || output <= 0
disp("Illegal value. Call with 0 < output < 1");
return
end
location = input .* ((1 - output) ./output).^ (1./exponent);
end
% The symbol for the subjective-price function in the flow diagrams is a
% hockey-stick-shaped curve in a box.
%
% Solomon, R. B., Conover, K., & Shizgal, P. (2017).
% Valuation of opportunity costs by rats working for rewarding electrical brain stimulation.
% PLOS ONE, 12(8), e0182120. https://doi.org/10.1371/journal.pone.0182120
function Psub = PsubFun(Pobj, PsubBend, PsubMin) % Subjective-price function
Psub = PsubMin + PsubBend .* (log(1 + exp((Pobj - PsubMin)./PsubBend)));
end
% The equation for PsubBsFun will return -Inf when Psub = PsubMin and a complex number when Psub includes one more values less
% than PsubMin. When Psub includes values between PsubMin and Pobj_0 (the value of Pobj at which Psub = 0),
% PsubBsFun returns negative values. The definition of PsubBsFun has thus been couched so as to return real values
% > 0 when Psub > Pobj_0, 0 when Psub = Pobj_0, and NaN otherwise. The reason that 0 is returned when obj_0, 0 when Psub = Pobj_0
% is to ensure consistency between the forward and backward solutions.
% Last revised by Peter Shizgal on 2019-09-08 15:06
function Pobj = PsubBsFun(Psub, PsubBend, PsubMin) % Back-solution of subjective-price function
Psub_Pobj_0 = PsubMin + PsubBend .* log(1 + exp((-PsubMin)./PsubBend));
above = Psub > Psub_Pobj_0;
eq = Psub == Psub_Pobj_0;
below = Psub < Psub_Pobj_0;
Pobj = PsubMin + PsubBend .* log(-1 + exp((Psub(above) - PsubMin)./PsubBend));
% Output will be real because Psub(above) > the value that drives Pobj to zero
Pobj(eq) = 0;
Pobj(below) = NaN;
end
% The next four functions are defined for consistency in formatting when composing more complex functions from the
% five building blocks (FilterFun, PsubFun, LogistNormFun, ScalarDivFun, ScalarMultFun) and their
% six back-solutions (FilterFunBS, PsubBsFun, LogistNormBsFun, LogistNormBsLocFun, ScalarDivBsFun, ScalarMultBsFun).
% Obviously, ScalarDiv function == "/", ScalarDivBsFun == "*", ScalarMultFun == "*", and
% ScalarMultBsFun == "/".
function quotient = ScalarDivFun(dividend, divisor)
quotient = (dividend ./ divisor);
end
function dividend = ScalarDivBsFun(divisor, quotient)
% Back solution of the scalar-division function (scalar multiplication)
dividend = (quotient .* divisor);
end
function product = ScalarMultFun(multiplicand, multiplier)
product = (multiplicand .* multiplier);
end
% Scalar multiplication is represented in the flow diagrams in three ways:
% 1) as a product operator (Pi notation) in a box (to compute the aggregate firing rate),
% 2) as a right-facing triangle (representing amplification), and
% 3) as a standard multiplication sign
function multiplicand = ScalarMultBsFun(multiplier, product)
% Back solution of the scalar-multiplication function (scalar division)
multiplicand = (product ./ multiplier);
end
Functions composing the reward-mountain model
% Dummy function to relate the duration of the burst of increased firing in the directly stimulated neurons to
% the pulse frequency. We set the two to be equal, a reasonable assumption if frequency-following fidelity is
% high.
function Dburst = fD(Dtrain)
Dburst = Dtrain;
end
% Function to compute the aggregate rate of firing required to produce a reward of half-maximal intensity
% see: Sonnenschein, B., Conover, K., & Shizgal, P. (2003).
% Growth of brain stimulation reward as a function of duration and stimulation strength.
% Behavioral Neuroscience, 117(5), 978–994. http://psycnet.apa.org/journals/bne/117/5/978
function FFaggHM = fH(C, D, RhoPi)
FFaggHM = ScalarMultFun(...
RhoPi,...
(1 + ScalarDivFun(...
C,...
D...
)...
)...
);
end
% Function to compute the average rate of firing required to produce a reward of half-maximal intensity
% see: Sonnenschein, B., Conover, K., & Shizgal, P. (2003).
% Growth of brain stimulation reward as a function of duration and stimulation strength.
% Behavioral Neuroscience, 117(5), 978–994. http://psycnet.apa.org/journals/bne/117/5/978
function FFhm = fHfiring(C, D, N, RhoPi, varargin)
if size(varargin,2) == 0
FFhm = ScalarDivFun(...
fH(...
C, ...
D, ...
RhoPi...
),...
N...
);
else
Kin = varargin{1};
FFhm = ScalarDivFun(...
fH(...
C, ...
D, ...
ScalarDivFun(...
RhoPi,...
Kin...
)...
),...
N...
);
end
end
% Function to compute the average pulse frequency required to produce a reward of half-maximal intensity
% see: Sonnenschein, B., Conover, K., & Shizgal, P. (2003).
% Growth of brain stimulation reward as a function of duration and stimulation strength.
% Behavioral Neuroscience, 117(5), 978–994. http://psycnet.apa.org/journals/bne/117/5/978
function FpulseHM = FpulseHMfun(C, D, Fbend, FPmax, Fro, N, RhoPi, varargin)
if size(varargin,2) == 0
FpulseHM = FilterFunBS(...
Fbend, ...
Fro, ...
fHfiring(C, ...
D, ...
N, ...
RhoPi),...
FPmax...
);
else
Kin = varargin{1};
FpulseHM = FilterFunBS(...
Fbend, ...
Fro, ...
fHfiring(C, ...
D, ...
N, ...
RhoPi,...
Kin...
),...
FPmax...
);
end
end
% Function to find the pulse frequency that produces a reward of half-maximal intensity given power growth of
% reward intensity
% FFmax, the maximal induced firing rate is obtained by setting the pulse frequency to 1000
function FPhmPG = find_FhmPG(FPmax, Fbend, Fro, g, Kin, Kout, RmaxPG)
FFmax = FilterFun(FPmax,Fbend,Fro);
FFhmPG = Kin .* ((0.5 .* RmaxPG) ./ Kout).^(1./g);
FPhmPG = FilterFunBS(Fbend, Fro, FFhmPG, FPmax);
end
% Function to translate the effect of an electrical input to the directly stimulated mfb stage
% of the series circuit model into an equivalent optical input to the dopaminergic stage.
% Both inputs produce the same level of firing in the dopamine neurons.
function FoptEquiv = FoptEquivFun(Felec, FmfbBend, FmfbHM, FmfbRO, gMFB, KrgE)
FoptEquiv = fRbsr(...
Felec,...
FmfbBend,...
FmfbHM,...
FmfbRO,...
gMFB,...
KrgE);
end
% Function to compute the equivalent optical drive on midbrain dopamine neurons
% due to trans-synaptic input from electrically excited MFB neurons
function FmfbDAdrive = FmfbDAdriveFun(FmfbDAbend, Felec, FFmfbDAmax, FoptEquiv, FmfbDAro)
FmfbDAdrive = ScalarMultFun(...
ScalarDivFun(...
FilterFun(...
Felec,...
FmfbDAbend,...
FmfbDAro...
),... % Filtered firing rate in MFB input
FFmfbDAmax...
),... % Proportion of maximum firing rate achieved
FoptEquiv...
); % Maximum equivalent optical pulse frequency applied to the dopamine neurons
end
% Find the electrical pulse frequency that produces half-maximal, summated reward intensity
% The solution is obtained by means of interpolation
%function [FsumHM, RsumMax] = find_FhmSum(Felec, Rsum)
function [FsumHM, RsumMax, CindAbove, CindBelow] = find_FhmSum(Felec, Rsum)
RsumMax = max(Rsum, [], 2); % max along each row
Rcrit = RsumMax./2;
RsumHMlog = Rsum > Rcrit;
RsumHMlogDiff = diff(RsumHMlog==1,1,2); % indices (offset by 1) marking transition from zero to one
% RsumHMlogDiff = diff(RsumHMlog); % indices (offset by 1) marking transition from zero to one
siz = size(RsumHMlogDiff);
IndDiffAbove = find(RsumHMlogDiff==1); % linear indices for first elements > Rcrit
[Rind Cind] = ind2sub(siz,IndDiffAbove);
% row & column indices for first elements > Rcrit offset by -1
CindSort(Rind) = Cind; % column indices sorted by row number of Rsum matrix
CindAbove = CindSort + 1; % row index of value just above Fh
CindBelow = CindSort;
% find exact Fhm by interpolation
siz = size(Rsum); % one more column than diff matrix
IndAbove = sub2ind(siz,1:length(Rind),CindAbove);
IndBelow = IndAbove - siz(1);
% The indices increment down the columns, so adjacent row indices are separated by the number of rows
RsumAbove = Rsum(IndAbove);
RsumBelow = Rsum(IndBelow);
slp = (RsumAbove - RsumBelow) ./ (Felec(CindAbove) - Felec(CindBelow));
FsumHM = (((Rcrit' - RsumBelow) ./ slp) + Felec(1,CindBelow))';
% All rows of Felec are the same
% Rcrit is a column vector but RsumBelow and Felec are row vectors
% Return FsumHM as a column vector so that it has the same dimensions as RsumMax
end
% Function to compute the normalized value of the reward-growth function
% F is a pulse frequency, whereas FFhm is a firing frequency
function RbsrNorm = fRbsrNorm(F, Fbend, FFhm, Fro, g)
RbsrNorm = LogistNormFun(...
g,...
FilterFun(...
F,...
Fbend,...
Fro...
),...
FFhm...
);
end
% Function to compute the scaled value of the reward-growth function
% F is a pulse frequency, whereas Fhm is a firing frequency
function Rbsr = fRbsr(F, Fbend, Fhm, Fro, g, Krg)
Rbsr = ScalarMultFun(...
Krg,...
LogistNormFun(...
g,...
FilterFun(...
F,...
Fbend,...
Fro...
),...
Fhm...
)...
);
end
function Rbsr = fRbsrFull(C, D, F, Fbend, Fro, g, Krg, N, RhoPi, varargin)
if size(varargin,2) == 0
Rbsr = ScalarMultFun(...
Krg,...
LogistNormFun(...
g,...
FilterFun(...
F,...
Fbend,...
Fro...
),...
fHfiring(...
C, ...
D, ...
N, ...
RhoPi...
)...
)...
);
else
Kin = varargin{1};
Rbsr = ScalarMultFun(Krg,...
LogistNormFun(g,...
FilterFun(...
F,...
Fbend,...
Fro...
),...
fHfiring(...
C, ...
D, ...
N, ...
RhoPi,...
Kin...
)...
)...
);
end
end
% Back solution of the full reward-growth function for the series-circuit (sc) model
% KrgE translates the normalized output of the upstream reward-growth function into
% the equivalent of an optical pulse frequency applied to the dopamine neurons
function FscHM = FscHMbs(FhmDA, FmfbBend, FhmMFB, FmfbRO, FPmax, gMFB, KrgEq)
FscHM = FilterFunBS(... % Electrical pulse frequency that produces HM reward intensity
FmfbBend,...
FmfbRO,...
LogistNormBsFun(... % MFB firing-rate that produces HM reward intensity
gMFB,... % exponent of MFB RG function
FhmMFB,... % location parameter of MFB RG function
ScalarMultBsFun(... % normalized output of MFB RG function
KrgEq,...
FhmDA...
)... % normalized output of MFB RG function
),... % MFB firing-rate that produces HM reward intensity
FPmax... % pulse frequency sufficient to maximize firing
); % Electrical pulse frequency that produces HM reward intensity
end
% Power function to compute reward growth
% Kin scales the firings, and the product is raised to the power of g
function R = fRpg(F, Fbend, Fro, g, Kin, Kout)
R = ScalarMultFun(...
Kout,...
ScalarDivFun(...
FilterFun(F,...
Fbend,...
Fro...
),...
Kin...
)...
.^g...
);
end
% Power function to compute reward growth, normalized to output values between 0 & 1
% A pulse frequency of 1000 is used to estimate Rmax
function Rnorm = RGpgNormFun(F, Fbend, Fro, g, Kin, Kout)
Rnorm = ScalarDivFun(...
fRpg(F, Fbend, Fro, g, Kin, Kout),...
fRpg(1000, Fbend, Fro, g, Kin, Kout)...
);
end
% Dummy function to compute subjective probability (output = input)
% Breton, Y.-A., Conover, K., & Shizgal, P. (2014).
% The effect of probability discounting on reward seeking: a three-dimensional perspective.
% Frontiers in Behavioral Neuroscience, 8, 284. https://doi.org/10.3389/fnbeh.2014.00284
% found that subjective and objective reward probabilites in an eICSS experiment
% were indistiguishable when the objective probability was 0.5 or greater.
% Thus, subjective and objective probabilities are equated here.
% We do not know at what point this equivalence breaks down, so use of
% this function is invalid for objective probabilities <0.5.
function ProbSub = ProbSubFun(ProbObj)
ProbSub = ProbObj;
end
% Function to compute subjective effort cost (output = input)
% This function is certainly not scalar! However, we haven't measured it yet.
% We make the assumption that the rate of subjective exertion is held constant
% under the conditions of our experiment, and we thus impose a fixed value on the subjective rate of exertion.
% Due to factors such as fatigue, this can't be entirely right. We hope
% that violations of this assumption aren't serious.
% To accommodate the hypothesis that subjective effort costs are modulated by dopamine tone,
% the optional KeffMod parameter scales the subjective effort cost.
% We use this function to set the subjective rates of exertion for both work and leisure.
function dotPhiSubScaled = fphiSc(dotPhiSub, Kec, varargin)
if size(varargin,2) == 0
dotPhiSubScaled = ScalarMultFun(...
dotPhiSub,...
Kec...
);
else
KeffMod = varargin{1};
dotPhiSubScaled = ScalarMultFun(...
dotPhiSub,...
ScalarDivFun(...
Kec,...
KeffMod...
)...
);
end
end
% Dummy function to compute subjective value of alternate ("leisure") activities (output = input)
% Again, we impose a fixed value, faute de mieux.
function RaaScaled = fRaa(Raa,Kaa)
RaaScaled = ScalarMultFun(...
Raa,...
Kaa...
);
end
% Function to compute the subjective equivalent of the objective-price location parameter
function PsubE = PsubEfun(dotPhiObj, Kaa, Kec, Krg, Raa, pObj, RnormMax, varargin)
if size(varargin,2) == 0
PsubE = ScalarDivFun(...
ScalarMultFun(...
ProbSubFun(pObj),...
ScalarMultFun(... % compute Rmax
Krg,...
RnormMax...
)...
),...
ScalarMultFun(...
fphiSc(dotPhiObj, Kec),...
fRaa(Raa, Kaa)...
)...
);
else
KeffMod = varargin{1};
PsubE = ScalarDivFun(...
ScalarMultFun(...
ProbSubFun(pObj),...
ScalarMultFun(... % compute Rmax
Krg,...
RnormMax...
)...
),...
ScalarMultFun(...
fphiSc(dotPhiObj, Kec, KeffMod),...
fRaa(Raa, Kaa)...
)...
);
end
end
% Function to compute the subjective-price location parameter for the power-growth model
function PsubEpg = PsubEpgFun(dotPhiObj, Kaa, Kec, Krg, ObjAA, pObj, Rmax)
PsubEpg = ScalarDivFun(...
ScalarMultFun(...
Rmax,...
ProbSubFun(pObj)...
),...
ScalarMultFun(...
fphiSc(dotPhiObj, Kec),...
fRaa(ObjAA, Kaa)...
)...
);
end
% Function to compute the objective-price location parameter
% This is the back-solution of the function that computes PsubE
function PobjE = PobjEfun(dotPhiObj, Kaa, Kec, Krg, ObjAA, pObj, PsubBend, PsubMin, RnormMax, varargin)
if size(varargin,2) == 0
PobjE = PsubBsFun(...
PsubEfun(dotPhiObj, Kaa, Kec, Krg, ObjAA, pObj, RnormMax), ...
PsubBend, ...
PsubMin...
);
else
KeffMod = varargin{1};
PobjE = PsubBsFun(...
PsubEfun(dotPhiObj, Kaa, Kec, Krg, ObjAA, pObj, RnormMax, KeffMod), ...
PsubBend, ...
PsubMin...
);
end
end
% Function to compute time allocation
% This function builds the mountain surface
% The values of the location parameters, FpulseHM & PobjE, are inputs
% These must be computed before this function is called
function T = TAfun(a, F, Fbend, FpulseHM, Fro, g, Pobj, PobjE, PsubBend, PsubMin, RnormMax)
T = LogistNormFun(...
a, ...
fRbsrNorm(...
F, ...
Fbend, ...
FpulseHM, ...
Fro, ...
g...
), ...
ScalarDivFun(...
ScalarMultFun(...
PsubFun(Pobj, ...
PsubBend, ...
PsubMin),...
RnormMax...
),...
PsubFun(PobjE, ...
PsubBend, ...
PsubMin)...
)...
);
end
% Function to compute time allocation using power growth of reward intensity
% This function builds the mountain surface
% The values of the location parameter, PsubEpg, is an input and must be computed before this function is called
function T = TApgFun(a, F, Fbend, Fro, g, Kin, Kout, Pobj, PsubBend, PsubEpg, PsubMin)
T = LogistNormFun(...
a, ...
RGpgNormFun(F, Fbend, Fro, g, Kin, Kout),...
ScalarDivFun(...
PsubFun(Pobj, ...
PsubBend, ...
PsubMin),...
PsubEpg...
)...
);
end
% Function to compute time allocation for the convergence model
% The reward-growth function must be computed beforehand. Its output, RGvec, must be a column vector.
function T = TAsumFun(a, Pobj, PobjE, PsubBend, PsubMin, RGvec, RnormMax)
T = LogistNormFun(...
a, ...
RGvec, ...
ScalarDivFun(...
ScalarMultFun(...
PsubFun(Pobj, ...
PsubBend, ...
PsubMin),...
RnormMax...
),...
PsubFun(PobjE, ...
PsubBend, ...
PsubMin)...
)...
);
end
Functions that plot graphs and set attributes
% Plot of the frequency-following function
% function FF_graph = plot_freqFoll(F, FR, fnroot, graphs2files, figdir, varargin)
function FF_graph = plot_freqFoll(F, FR, fnroot, varargin)
global graphs2files FigDir
if size(varargin,2) == 0
% To accommodate multiple plots based on input matrices, must compute {min, max} in 2 stages
% Revise this section to operate in the log domain and then transform back to linear
xmin = min(min(F)) - (0.1*(max(max(F))-min(min(F))));
xmax = max(max(F)) + (0.1*(max(max(F))-min(min(F))));
ymin = min(min(FR)) - (0.1*(max(max(FR))-min(min(FR))));
ymax = max(max(FR)) + (0.1*(max(max(FR))-min(min(FR))));
else
xmin = varargin{1};
xmax = varargin{2};
ymin = varargin{3};
ymax = varargin{4};
end
FF_graph = figure;
FF_graph.Visible = 'off';
FF_graph.Position = [0 0 600 600];
gmfbfr = loglog(F,FR);
% gmfbfr.LineWidth = 4;
arrayfun(@(x) set(x,'LineWidth',4), gmfbfr); % Needed if plot includes multiple lines
grid on;
ax1 = gca;
ax1.XLim = [xmin xmax];
ax1.YLim = [ymin ymax];
ax1.XLabel.String = 'Pulse Frequency_{ }';
ax1.YLabel.String = 'Firing Frequency_{ }';
ax1.FontSize = 24;
ax1.XLabel.FontAngle = 'italic';
ax1.YLabel.FontAngle = 'italic';
% There is an unsolved problem in retrieval of tick values. This works when stepping through the code but
% not when running it! The following is a work-around for setting tick and tick-label values.
min_xtick = ceil(log10(xmin)); % power of ten above minimum x value
max_xtick = floor(log10(xmax)); % power of ten below maximum x value
log_xt = min_xtick:max_xtick;
xt = 10.^log_xt;
num_xtick = length(xt);
xtick_lab = cell(num_xtick,1);
for j = 1:num_xtick
xtick_lab{j} = sprintf('%4.0f',xt(j));
end
xticks(xt);
xticklabels(xtick_lab);
min_ytick = ceil(log10(ymin)); % power of ten above minimum x value
max_ytick = floor(log10(ymax)); % power of ten below maximum x value
log_yt = min_ytick:max_ytick;
yt = 10.^log_yt;
num_ytick = length(yt);
ytick_lab = cell(num_ytick,1);
for j = 1:num_ytick
ytick_lab{j} = sprintf('%4.0f',yt(j));
end
yticks(yt);
yticklabels(ytick_lab);
if graphs2files
saveas(gca,fullfile(FigDir,strcat(fnroot,'_freq_resp.fig')));
saveas(gca,fullfile(FigDir,strcat(fnroot,'_freq_resp.png')));
end
end
% Plot of the strength-duration function for trains
function SDT_graph = plot_StrDurTrains(Delec,FelecHM)
global graphs2files FigDir
SDT_graph = figure;
SDT_graph.Visible = 'off';
SDT_graph.Position = [0 0 600 600];
gmfbfr = loglog(Delec,FelecHM);
gmfbfr.LineWidth = 4;
grid on;
ax1 = gca;
ymax = max(FelecHM)*10^0.1;
ymin = min(FelecHM)*10^-0.1;
ax1.YLim = [ymin ymax];
ax1.XLabel.String = 'Train Duration (s)_{ }';
ax1.YLabel.String = 'F_{hm }';
ax1.FontSize = 24;
ax1.XLabel.FontAngle = 'italic';
ax1.YLabel.FontAngle = 'italic';
ax1.XTickLabels = {'0.1', '1', '10', '100'};
if graphs2files
saveas(gca,fullfile(FigDir,'mfb_StrDurTrains.fig'));
saveas(gca,fullfile(FigDir,'mfb_StrDurTrains.png'));
end
end
% Plot one or more reward-growth functions
function out_graph = plot_RG(Fmat,RMat,pnam,pVec,fnam,TitleStr,linlog, varargin)
global graphs2files FigDir
fh = figure;
fh.Position = [0 0 600 600];
fh.Visible = 'off';
plab = regexprep(regexprep(pnam,'\W',''),'_','');
switch linlog
case 'lin' % linear y-axis, logarithmic x-axis
gRG = semilogx(Fmat,RMat);
ymin = -0.05;
ymax = max(max(RMat)) * 1.05;
if graphs2files
if size(varargin,2) == 0
% figfile = fullfile(FigDir,char(strcat('RG_',fnam','_','semilog.fig')));
% pngfile = fullfile(FigDir,char(strcat('RG_',fnam','_','semilog.png')));
figfile = fullfile(FigDir,char(strcat('RG_',fnam,'_','semilog.fig')));
pngfile = fullfile(FigDir,char(strcat('RG_',fnam,'_','semilog.png')));
else
name_suffix = varargin{1};
figfile = ...
fullfile(FigDir,char(strcat('RG_',fnam','_',name_suffix,'_','semilog.fig')));
pngfile = ...
fullfile(FigDir,char(strcat('RG_',fnam','_',name_suffix,'_','semilog.png')));
end
end
case 'log' % logarithmic y-axis, logarithmic x-axis
gRG = loglog(Fmat,RMat);
ymin = 0.01;
ymax = max(max(RMat)) * 10^0.1;
if graphs2files
if size(varargin,2) == 0
figfile = fullfile(FigDir,char(strcat('RG_',fnam','_','loglog.fig')));
pngfile = fullfile(FigDir,char(strcat('RG_',fnam','_','loglog.png')));
else
name_suffix = varargin{1};
figfile = ...
fullfile(FigDir,char(strcat('RG_',fnam','_',name_suffix,'_','loglog.fig')));
pngfile = ...
fullfile(FigDir,char(strcat('RG_',fnam','_',name_suffix,'_','loglog.png')));
end
end
end
for j = 1:length(gRG)
gRG(j).LineWidth = 4;
end
grid on;
ax1 = gca;
ax1.Position = [0.175 0.15 0.75 0.75];
ax1.XLim = [10 1000];
ax1.YLim = [ymin ymax];
ax1.XLabel.String = 'Pulse Frequency_{ }';
ax1.YLabel.String = 'Reward Intensity_{ }';
ax1.FontSize = 24;
ax1.XLabel.FontAngle = 'italic';
ax1.YLabel.FontAngle = 'italic';
ax1.XTick = [1 10 100 1000];
ax1.XTickLabels = {'1', '10', '100', '1000'};
ax1.LineWidth = 1;
switch linlog
case 'lin'
XlabPos = ax1.XLabel.Position;
XlabPos(2) = ax1.YLim(1) - ((ax1.YLim(2) - ax1.YLim(1)) * 0.075); % 7.5% below lower limit
ax1.XLabel.Position = XlabPos;
case 'log'
XlabPos = ax1.XLabel.Position;
XlabPos(2) = 10^(log10(ax1.YLim(1)) - ((log10(ax1.YLim(2)) - log10(ax1.YLim(1)))*0.075)); % 7.5% below lower limit
ax1.XLabel.Position = XlabPos;
end
YlabPos = ax1.YLabel.Position;
YlabPos(1) = 10^(log10(ax1.XLim(1)) - ((log10(ax1.XLim(2)) - log10(ax1.XLim(1)))*0.1)); % 10% below lower limit
ax1.YLabel.Position = YlabPos;
% if min(pVec) < 0.01
% formatstring = '%3.2e';
% else
% formatstring = '%0.2f';
% end
if min(pVec) < 10
formatstring = '%3.2f';
else
formatstring = '%3.0f';
end
% Make a cell array of the legend values and pass to the legend function
lgnd = legend(arrayfun(@(x) sprintf(formatstring,x), pVec, 'UniformOutput', false),...
'Location', 'best');
lgnd.Title.String = strjoin(pnam,'\n');
title(ax1, TitleStr, 'FontSize', 30);
out_graph = gcf;
if graphs2files
saveas(gca,figfile);
saveas(gca,pngfile);
end
end
% Plot of the subjective-price function
function SP_graph = plot_PsubFun(Pobj,Psub,graphs2files, figdir)
SP_graph = figure;
SP_graph.Visible = 'off';
SP_graph.Position = [0 0 600 600];
gPsub = loglog(Pobj,Psub);
gPsub.LineWidth = 4;
grid on;
ax1 = gca;
ax1.XLim = [0.1 100];
ax1.YLim = [0.5 100];
ax1.XTickLabels = {'0.1', '1', '10', '100'};
ax1.YTickLabels = {'1', '10', '100'};
ax1.LineWidth = 1;
ax1.XLabel.String = 'Objective Price (s)_{ }';
ax1.YLabel.String = 'Subjective Price (s)_{ }';
ax1.FontSize = 24;
ax1.XLabel.FontAngle = 'italic';
ax1.YLabel.FontAngle = 'italic';
if graphs2files
saveas(gca,fullfile(figdir,'Subj_price.fig'));
saveas(gca,fullfile(figdir,'Subj_price.png'));
end
end
% Function to plot a single mountain
% F is a logarithmically spaced vector of pulse-frequencies
% Pobj is a logarithmically spaced vector of prices
% T is a square matrix of time-allocation values computed with price as the row variable
% and pulse frequency as the column variable.
% Visible is a character vector ('on', or 'off') that determines whether or not the figure is displayed.
% The optional arguments are xmin, xmax, ymin, ymax.
function gh3d2 = plot_MTN(F, Pobj, T, Visible, mtn_root, title_str, ...
graphs2files, figdir, varargin)
logF = log10(F);
logP = log10(Pobj);
[X,Y] = meshgrid(logP, logF);
Z = T;
if size(varargin,2) == 0
xmin = 0;
xmax = 2.6;
ymin = 1;
ymax = 3;
else
xmin = varargin{1};
xmax = varargin{2};
ymin = varargin{3};
ymax = varargin{4};
end
gh3d = figure;
gh3d.Visible = Visible;
gh3d.Position = [0 0 600 600];
surface(X, Y, Z, 'CData', Z,...
'FaceLighting','gouraud',...
'FaceColor',[0.5 0.5 0.5], 'FaceAlpha', 0.5,...
'EdgeColor', [0.5 0.5 0.5],...
'SpecularColorReflectance',0.5, 'SpecularExponent',4, 'SpecularStrength',1);
light('Parent',gca,'Position',[10 -5 2]);
view([40 30]);
ax3d = findall(gca,'Type','Axes');
title(ax3d, title_str, 'FontSize', 36);
ax3d.Position = [0.15, 0.15, 0.75, 0.75];
ax3d.FontSize = 16;
ax3d.XLabel.String = 'Objective Price (s)';
ax3d.XLabel.FontSize = 22;
ax3d.XLabel.Rotation = -27;
ax3d.XLabel.Units = 'Normalized';
ax3d.XLabel.Position = [0.325 -0.0625 0];
ax3d.XAxis.Scale = 'linear';
ax3d.XLim = [xmin,xmax];
ax3d.XTick = [0, 0.5, 1, 1.5, 2, 2.5];
ax3d.XTickLabel = {'1.0','3.2','10','32','100','316'};
ax3d.YLabel.String = 'Pulse frequency (pps)';
ax3d.YLabel.FontSize = 22;
ax3d.YLabel.Rotation = 32;
ax3d.YLabel.Units = 'Normalized';
ax3d.YLabel.Position = [0.85 0.05 0];
ax3d.YAxis.Scale = 'linear';
ax3d.YLim = [ymin,ymax];
ax3d.YTick = [0, 0.5, 1, 1.5, 2, 2.5, 3];
ax3d.YTickLabel = {'1','3.2','10','32', '100', '316', '1000'};
ax3d.ZLim = [0,1.1];
ax3d.ZLabel.String = 'Time Allocation';
ax3d.ZLabel.FontSize = 22;
ax3d.ZLabel.Rotation = 91;
ax3d.ZLabel.Units = 'Normalized';
ax3d.ZLabel.Position = [-0.125 0.475 0];
ax3d.DataAspectRatio = [1 0.75 0.8];
hold on;
z_shim = 0.01;
z_shim = Z + z_shim;
contvec = 0.1:0.1:0.9;
[~, ch] = contour3( X, Y, z_shim, contvec);
colormap('jet');
set(ch, 'LineWidth', 2);
gh3d2 = figure(gcf);
gh3d2.Visible = Visible;
grid on;
hold off;
if graphs2files
saveas(gca,fullfile(figdir,strcat(mtn_root,'.fig')));
saveas(gca,fullfile(figdir,strcat(mtn_root,'.png')));
end
end
% Function to plot a single contour graph
% F is a logarithmically spaced vector of pulse-frequencies.
% Pobj is a logarithmically spaced vector of prices.
% T is square matrix of time-allocation values computed with price as the row variable
% and pulse frequency as the column variable.
% Visible is a character vector ('on', or 'off') that determines whether or not the figure is displayed.
% The optional arguments are xmin, xmax, ymin, ymax.
function cont1 = plot_contour(F, Pobj, T, Pobj_e, Fhm, Visible, mtn_root, ...
title_str, annot_str, graphs2files, figdir, varargin)
logF = log10(F);
logP = log10(Pobj);
[X,Y] = meshgrid(logP, logF);
Z = T;
if size(varargin,2) == 0
xmin = 0;
xmax = 2.6;
ymin = 1;
ymax = 3;
else
xmin = varargin{1};
xmax = varargin{2};
ymin = varargin{3};
ymax = varargin{4};
end
contvec = 0.1:0.1:0.9;
cont1 = figure;
cont1.Visible = Visible;
cont1.Position = [0 0 600 600];
contour(X,Y,Z,contvec,'LineWidth',1.25);
colormap('jet');
hold on;
axc1 = findall(gca,'Type','Axes');
title(axc1, title_str, 'FontSize', 36);
axc1.Position = [0.15, 0.175, 0.825, 0.75];
axc1.FontSize = 18;
axc1.XLim = [xmin,xmax];
axc1.XLabel.String = 'Objective Price (s)';
axc1.XLabel.FontSize = 24;
% Adjust position of X-axis label
xlpos = axc1.XLabel.Position;
xlpos(2) = ymin - ((ymax - ymin) * 0.1);
axc1.XLabel.Position = xlpos;
axc1.XTick = [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5];
axc1.XTickLabel = {'','','0.1','0.32','1.0','3.2','10','32','100', '316'};
axc1.YLim = [ymin,ymax];
axc1.YLabel.String = 'Pulse Frequency (pps)';
axc1.YLabel.FontSize = 24;
% Adjust position of Y-axis label
ylpos = axc1.YLabel.Position;
ylpos(1) = xmin - ((xmax - xmin) * 0.125);
axc1.YLabel.Position = ylpos;
axc1.YTick = [0. 0.5, 1, 1.5, 2, 2.5, 3];
axc1.YTickLabel = {'1', '3.2', '10', '31.6', '100', '316', '1000'};
FhmLine1 = line('XData', axc1.XLim,'YData', [log10(Fhm), log10(Fhm)]);
FhmLine1.Color = 'r';
FhmLine1.LineWidth = 4;
FhmLine1.LineStyle = '--';
PeLine1 = line('XData', [log10(Pobj_e), log10(Pobj_e)],'YData', axc1.YLim);
PeLine1.Color = 'b';
PeLine1.LineWidth = 4;
PeLine1.LineStyle = '--';
tb2 = annotation('textbox');
tb2.String = annot_str;
tb2.FontSize = 18;
tb2.Position = [0.7, 0.25, 0.15, 0.075];
hold off;
legend([FhmLine1, PeLine1], [{'F_{hm }'}, {'P_e '}], 'Location','northwest');
if graphs2files
saveas(gca,fullfile(figdir,strcat(mtn_root,'.fig')));
saveas(gca,fullfile(figdir,strcat(mtn_root,'.png')));
end
end
% Function to produce a bargraph showing changes in the location parameters
function bg = plot_bg(Fhm1, Fhm2, Pobj_e1, Pobj_e2, Visible, bg_root,...
graphs2files, figdir, varargin)
logFhmShift = log10(Fhm2) - log10(Fhm1);
logPobjShift = log10(Pobj_e2) - log10(Pobj_e1);
if size(varargin,2) == 0
span = max(logFhmShift, logPobjShift) - min(logFhmShift, logPobjShift);
ymax = max(logFhmShift, logPobjShift) + 0.2 * span;
ymin = min(logFhmShift, logPobjShift) - 0.2 * span;
else
ymin = varargin{1};
ymax = varargin{2};
end
bg = figure;
bg.Visible = Visible;
bg.Position = [0 0 600 800];
bar(1,logFhmShift, 'r');
hold on;
bar(2,logPobjShift, 'b');
hold off;
ax = findall(gca,'Type','Axes');
ax.Position = [0.2 0.15 0.6 0.8]; %[left bottom width height]
ax.XLim = [0.5,2.5];
ax.YLim = [ymin, ymax];
ax.YLabel.String = 'Shift';
ylabpos = ax.YLabel.Position;
ylabpos(1) = 0.15;
ax.YLabel.Position = ylabpos;
ax.XAxis.FontSize = 22;
ax.YAxis.FontSize = 22;
ax.YGrid = 'on';
ax.XTick = [1,2];
ax.XTickLabel = {'log_{10 }(F_{hm })', 'log_{10 }(P_e )'};
bh = findall(gcf,'Type','Bar');
bh(1).BarWidth = 0.6;
bh(2).BarWidth = 0.6;
if graphs2files
saveas(gca,fullfile(figdir,strcat(bg_root,'.fig')));
saveas(gca,fullfile(figdir,strcat(bg_root,'.png')));
end
end
% Function to produce a bargraph showing changes in the location parameters
function bg = plot_bgStar(Fhm1, Fhm2, FhmStar1, FhmStar2, Pobj_e1, Pobj_e2, PobjEstar1, PobjEstar2, ...
Visible, bg_root, graphs2files, figdir, varargin)
logFhmShift = log10(Fhm2) - log10(Fhm1);
logPobjShift = log10(Pobj_e2) - log10(Pobj_e1);
logFhmStarShift = log10(FhmStar2) - log10(FhmStar1); % These are the corrected values
logPobjStarShift = log10(PobjEstar2) - log10(PobjEstar1); % These are the corrected values
if size(varargin,2) == 0
span = max(logFhmShift, logPobjShift) - min(logFhmShift, logPobjShift);
ymax = max(logFhmShift, logPobjShift) + 0.2 * span;
ymin = min(logFhmShift, logPobjShift) - 0.2 * span;
else
ymin = varargin{1};
ymax = varargin{2};
end
bwidth = 0.6; % width of bars
FhmStarX = [1-(bwidth/2), 1+(bwidth/2)];
FhmStarY = [logFhmStarShift, logFhmStarShift];
PeStarX = [2-(bwidth/2), 2+(bwidth/2)];
PeStarY = [logPobjStarShift, logPobjStarShift];
bg = figure;
bg.Visible = Visible;
bg.Position = [0 0 600 800];
bar(1,logFhmShift, 'r');
hold on;
bar(2,logPobjShift, 'b');
bh = findall(gcf,'Type','Bar');
bh(1).BarWidth = bwidth;
bh(2).BarWidth = bwidth;
% Add lines showing corrected location parameters
plot(FhmStarX, FhmStarY, 'Color','c', 'LineStyle', '-.', 'LineWidth',6);
plot(PeStarX, PeStarY, 'Color','c', 'LineStyle', '-.', 'LineWidth',6);
hold off;
ax = findall(gca,'Type','Axes');
ax.Position = [0.2 0.15 0.6 0.8]; %[left bottom width height]
ax.XLim = [0.5,2.5];
ax.YLim = [ymin, ymax];
ax.YLabel.String = 'Shift';
ylabpos = ax.YLabel.Position;
ylabpos(1) = 0.15;
ax.YLabel.Position = ylabpos;
ax.XAxis.FontSize = 22;
ax.YAxis.FontSize = 22;
ax.YGrid = 'on';
ax.XTick = [1,2];
ax.XTickLabel = {'log_{10 }(F_{hm })', 'log_{10 }(P_e )'};
bh = findall(gcf,'Type','Bar');
bh(1).BarWidth = 0.6;
bh(2).BarWidth = 0.6;
if graphs2files
saveas(gca,fullfile(figdir,strcat(bg_root,'.fig')));
saveas(gca,fullfile(figdir,strcat(bg_root,'.png')));
end
end
% Modify attributes of a 2D graph. handle_attribute and attribute_value can be cell arrays,
% in which case attribute value must have the same number of rows as the length of the handle
% and the same number of columns as the number on handle_attributes.
% The restrictions of the set command apply. For example, this function cannot modify legend attributes.
function out_graph = modify_2D_graph(in_graph, handle_type, handle_attribute, attribute_value, ...
gFileNam, graphs2files, figdir)
gh = findall(in_graph, 'Type', handle_type);
set(gh, handle_attribute, attribute_value);
if strcmp(handle_attribute,'XLim') || strcmp(handle_attribute,'YLim')
ax = findall(gcf, 'Type', 'Axes');
XlabPos = ax.XLabel.Position;
YlabPos = ax.YLabel.Position;
LinLogX = ax.XAxis.Scale;
switch LinLogX
case 'linear'
XlabPos(1) = ax.XLim(1) + ((ax.XLim(2) - ax.XLim(1)) / 2);
YlabPos(1) = ax.Xlim(1) - ((ax.XLim(2) - ax.XLim(1)) * 0.1); % 10% below lower limit
case 'log'
XlabPos(1) = 10^(log10(ax.XLim(1)) + ((log10(ax.XLim(2)) - log10(ax.XLim(1))) / 2));
YlabPos(1) = 10^(log10(ax.XLim(1)) - ((log10(ax.XLim(2)) - log10(ax.XLim(1)))*0.1)); % 10% below lower limit
end
ax.XLabel.Position = XlabPos;
LinLogY = ax.YAxis.Scale;
switch LinLogY
case 'linear'
XlabPos(2) = ax.YLim(1) - ((ax.YLim(2) - ax.YLim(1)) * 0.0/075); % 7.5% below lower limit;
YlabPos(2) = ax.YLim(1) + (ax.YLim(1) + ((ax.YLim(2) - ax.YLim(1)) / 2));
case 'log'
XlabPos(2) = 10^(log10(ax.YLim(1)) - ((log10(ax.YLim(2)) - log10(ax.YLim(1)))*0.075)); % 7.5% below lower limit
YlabPos(2) = 10^(log10(ax.YLim(1)) + ((log10(ax.YLim(2)) - log10(ax.YLim(1))) / 2));
end
ax.YLabel.Position = YlabPos;
end
out_graph = gcf;
if graphs2files
figfile = fullfile(figdir, strcat(gFileNam, '.fig'));
saveas(gcf,figfile);
pngfile = regexprep(figfile, '.fig','.png');
saveas(gcf,pngfile);
end
end
Functions that display multi-panel Matlab figures
function dual_sub = dual_subplot(g1, g2, dual_sub_out, graphs2files, figdir)
dual_sub = figure;
dual_sub.Visible = 'off';
dual_sub.Units = 'pixels';
dual_sub.Position = [0 0 1500 600];
colormap 'jet';
ax(1) = findall(g1, 'Type', 'Axes');
ax(2) = findall(g2, 'Type', 'Axes');
% Collect legend properties if a legend exists in the input figures
for j = 1:2
lg(j) = false;
if ~isempty(ax(j).Legend)
lg(j) = true;
lgh = ax(j).Legend;
lgstr(j).str = lgh.String;
lgstr(j).tstr = lgh.Title.String;
end
end
for j = 1:2
ax_copy(j) = copyobj(ax(j), dual_sub);
sph(j) = subplot(1,2,j,ax_copy(j));
% if ~isempty(ax(j).Legend)
if lg(j)
% splg(j) = legend(sph(j),lgstr(j).str,'Location','southeast');
% This position is rigid, but it should work in this specific case.
% The reward-growth graphs have legends; the mountain plots do not.
splg(j) = legend(sph(j),lgstr(j).str,'Location','best');
splg(j).Title.String = lgstr(j).tstr;
end
end
if graphs2files
saveas(dual_sub,fullfile(figdir,strcat(dual_sub_out,'.fig')));
saveas(dual_sub,fullfile(figdir,strcat(dual_sub_out,'.png')));
end
end
% function to plot the quad-panel display of the contour and bar graphs
% cont1 is plotted twice, once in the upper left and once in the lower right
function quad_sub = quad_subplot(cont1, cont2, bg, quad_sub_out, bg_root, ...
graphs2files, figdir)
quad_sub = figure;
quad_sub.Visible = 'off';
quad_sub.Units = 'pixels';
quad_sub.Position = [0 0 1500 1500];
colormap 'jet';
fh(1) = cont1;
fh(2) = bg;
fh(3) = cont2;
fh(4) = cont1;
for j = 1:4
ax(j) = findall(fh(j), 'Type', 'Axes');
end
% Collect textbox properties if a textbox exists in the input figure
for j = 1:4
tbh = findall(fh(j),'Type','TextBox');
tb(j) = false;
if ~isempty(tbh)
tb(j) = true;
tbstr(j) = tbh.String;
end
end
bg_str = findall(ax(2), '-Property', 'FontSize');
bg_str = findall(bg_str(:),'-Property','FontSize');
for j = 1:length(bg_str)
bg_str(j).FontSize = bg_str(j).FontSize * 1.75;
end
for j = 1:4
ax_copy(j) = copyobj(ax(j), quad_sub);
sph(j) = subplot(2,2,j,ax_copy(j));
if tb(j)
axpos = sph(j).Position;
x_offset = axpos(3) * 0.7;
y_offset = axpos(4) * -0.1;
tbdim = [(axpos(1) + x_offset) (axpos(2) + y_offset) 0.1 0.1];
ann(j) = annotation('textbox', tbdim, 'String', tbstr(j),'FitBoxToText','on');
ann(j).FontSize = 16;
end
end
if graphs2files
saveas(quad_sub,fullfile(figdir,strcat(quad_sub_out,'.fig')));
saveas(quad_sub,fullfile(figdir,strcat(quad_sub_out,'.png')));
end
end
function out_fig = adjust_right_panel(dual_sub, h, v)
out_fig = dual_sub;
out_fig.Visible = 'off';
ax = findall(out_fig,'Type','Axes');
ax1_pos = ax(1).Position;
ax1_ip = ax(1).InnerPosition;
ax1_op = ax(1).OuterPosition;
ax(1).InnerPosition = ax(1).OuterPosition; % tight borders
lgutter = ax1_op(1) - 0.5;
panel_width = 0.5 - lgutter;
panel_height = 1;
new_ax_width = ax1_pos(3) * h; % rescale height
new_ax_height = ax1_pos(4) * v; % rescale width
new_hmargin = (panel_width - new_ax_width) / 2;
new_vmargin = (panel_height - new_ax_height) / 2;
ax1_newxpos = 0.5 + lgutter + new_hmargin;
ax1_newpos(1)= ax1_newxpos;
ax1_newvpos = new_vmargin; % unnecessary, but included for consistency
ax1_newpos(2) = ax1_newvpos;
ax1_newpos(3) = new_ax_width;
ax1_newpos(4) = new_ax_height;
ax(1).Title.String = 'simulated data';
ax(1).Position = ax1_newpos;
end
Functions that create & display images from stored files
function show_imported_graphic(gnam, mag, ImpFigDir)
imshow(fullfile(ImpFigDir,gnam),'Border','tight','InitialMagnification',mag);
end
function F = make_fig_from_png(png_file, mag)
img = imread(png_file, 'png');
sz = size(img);
F = figure;
image(img);
pos = F.Position; % the conventional x and y dimensions are reversed
pos(3) = sz(2) * mag;
pos(4) = sz(1) * mag;
F.Position = pos;
axis tight;
ah = findall(F,'Type','Axes');
ah.Visible = 'off';
F.Visible = 'off';
end